Abstract. We study the uniform approximation of boundary layer functions exp(−x/d) for x ∈ (0, 1), d ∈ (0, 1], by the p and hp versions of the finite element method. For the p version (with fixed mesh), we prove super-exponential convergence in the range p + 1/2 > e/(2d). We also establish, for this version, an overall convergence rate of O(p −1 √ ln p) in the energy norm error which is uniform in d, and show that this rate is sharp (up to the √ ln p term) when robust estimates uniform in d ∈ (0, 1] are considered. For the p version with variable mesh (i.e., the hp version), we show that exponential convergence, uniform in d ∈ (0, 1], is achieved by taking the first element at the boundary layer to be of size O(pd).Numerical experiments for a model elliptic singular perturbation problem show good agreement with our convergence estimates, even when few degrees of freedom are used and when d is as small as, e.g., 10 −8 . They also illustrate the superiority of the hp approach over other methods, including a low-order h version with optimal "exponential" mesh refinement.The estimates established in this paper are also applicable in the context of corresponding spectral element methods.