This paper deals with the implementation of a MATLAB solver for singular boundary value problems with a singularity of the first kind. This solver includes a mesh adaption routine that selects meshes according to the local smoothness of the solution rather than to the smoothness of the direction field. The mesh selection is based on a new global error estimator for collocation schemes, which has recently been developed at the Institute for Applied Mathematics and Numerical Analysis and proves robust with respect to the singularity.