Abstract We consider the numerical solution of a singularly perturbed linear self-adjoint boundary value problem. Assuming that the coefficients of the differential equation are smooth, we construct and analyze finite difference methods that converge both with high order and uniformly with respect to the singular perturbation parameter. The analysis is done on a locally quasiuniform mesh, which permits its extension to the case of adaptive meshes which may be used to improve the solution. Numerical examples are presented to demonstrate the effectiveness of the method and its low computational cost. The convergence obtained in practice satisfies the theoretical predictions.