We describe a mesh selection strategy for the numerical solution of boundary value problems for singular ordinary differential equations. This mesh adaptation procedure is implemented in our MATLAB code \texttt{sbvp} which is based on polynomial collocation. We prove that under realistic assumptions our mesh selection strategy serves to approximately equidistribute the global error of the collocation solution, thus enabling to reach prescribed tolerances efficiently. Moreover, we demonstrate that this strategy yields a favorable performance of the code and compare its computational effort with other implementations of polynomial collocation.