We describe a robust, adaptive algorithm for the solution of singularly perturbed twopoint boundary value problems. Many different phenomena can arise in such problems, including boundary layers, dense oscillations, and complicated or ill-conditioned internal transition regions. Working with an integral equation reformulation of the original differential equation, we introduce a method for error analysis which can be used for mesh refinement even when the solution computed on the current mesh is underresolved. Based on this method, we have constructed a black-box code for stiff problems which automatically generates an adaptive mesh resolving all features of the solution. The solver is direct and of arbitrarily high-order accuracy and requires an amount of time proportional to the number of grid points.

Asymptotic and numerical methods are used to study several classes of singularly perturbed boundary value problems for which the underlying homogeneous operators have exponentially small eigenvalues. Examples considered include the familiar boundary layer resonance problems and some extensions, and certain linearized equations associated with metastable internal layer motion. For the boundary layer resonance problems, a systematic projection method, motivated by the work of De Groen [SIAM J. Math. Anal. 11, (1980), pp. 1-22], is used to analytically calculate high order asymptotic solutions. This method justifies and extends some previous results obtained from the variational method of Grasman and Matkowsky [SIAM J. Appl. Math. 32, (1977), pp. 588-597]. A numerical approach, based on an integral equation formulation, is used to accurately compute boundary layer resonance solutions and their associated exponentially small eigenvalues. For various examples, the numerical results are show...