We investigate mesh selection strategies based on global error estimation for the numerical solution of singular boundary value problems of the form

z'(t)=M(t)/t z(t)+f(t,z(t)), t∈ (0,1],	(1)
B0z(0)+B1z(1)=β,	(2)
z ∈ C[0,1].	(3)

where B₀, B₁ and β are suitably chosen to yield a well-posed problem. Collocation methods of variable order p are used to obtain the basic solution and different estimators of the global error are investigated. An estimator due to Zadunaisky works satisfactorily as in the classical case - for collocation schemes where no superconvergence effects are present - if only mesh points are used and collocation points neglected. However, this means that most of the computed information is discarded and that the step-sizes for the implicit Euler method which is involved in the error estimation are undesirably large. However, the same estimation procedure is not asymptotically correct if the collocation points are included in the process. Therefore, we propose a modification of the classical Zadunaisky method which makes it possible to obtain an asymptotically correct estimate in the collocation and mesh points. The information about the global error obtained thus is subsequently used for mesh selection. Our method turns out to yield meshes that are well suited to equidistribute the global error on the interval [0,1] and is not negatively affected by the singularity. This is an improvement upon classical mesh selection strategies which usually yield unnecessarily fine meshes near the singularity. The theory is supported by numerous numerical examples.