We revisit and extend the least-squares bubble function (LSBF) enrichment technique originally proposed by Yazdani and Nassehi, and derive and apply higher-order polynomials (degrees p = 3, 4) to a broader class of benchmark problems. The bubble function coefficients are determined by minimizing an element level L2 residual functional. An adaptive element-level order-selection rule based on the mesh Peclet and Damkohler numbers is proposed, drawing an analogy with local p-refinement in h-p finite element methods. The method is validated on three benchmark problems: (i) a singularly perturbed convection–diffusion–reaction (CDR) equation across four parameter regimes, (ii) a stiff two-point boundary value problem with a boundary layer of thickness O(ε) = 10e−4, and (iii) a Two-dimensional diffusion–reaction problem on a unit square. Systematic mesh-refinement studies confirm that LSBF at p = 3, 4 substantially outperforms standard Galerkin and is competitive with or superior to SUPG in reaction-dominated and mixed regimes. Some advantages and the limitations of the LSBF method are discussed.