Consider the singular ODE 


$y''+\frac{y'}{r}+p(r)y=0 \ \ with \ \ y(0)=1 \ \ and \ \ y'(0)=0$. 


Theoretically such solution exists and is unique if $p$ is nice. Is there a method to numerically solve this equation on some interval $[0, \alpha]$ with an explicit global error bound? I am not sure how to impose the initial condition at the singular point $r=0$.



Numerical solution of singular ODE
The [SLEIGN2](http://www.math.niu.edu/SL2/) software package can handle Sturm-Liouville problems with regular as well as singular end points. I believe you can find some papers on that web page that document various aspects of the methods used in the code.


