The numerical solution of singularly perturbed two-point boundary value problems in ordinary differential equations is considered. Implementation methods for general-purpose solvers of first-order linear systems are examined, with the basic difference scheme being collocation at Gaussian points. Adaptive mesh selection is based on localized error estimates at the collocation points These methods are implemented as modifications to the successful collocation code COLSYS, which was originally designed for mildly stiff problems only. Efficient high-order approximations to extremely stiff problems are obtained, and comparisons to COLSYS show that the modifications work relatively much better as the singular perturbation parameter gets small (i.e., the problem gets stiff), for both boundary-layer and turning-point problems.