(Why'd this get moved here, of all places? It's a math problem. Numerical methods, really.)

I guess this is a well-solved problem. Could anyone give me some help?

Well, more of an ill-conditioned problem, really. :)
I'm not at all a numerics expert, but the simplest (but surely not most efficient) numerical method that comes to mind would be to calculate the integral from h to infinity and shrink h until you're within convergence limits. You could also use Richardson extrapolation to improve on that.

The 1/r potential term is best handled by recognizing that the value from the analytical solution at r = 0 for the hydrogen atom is 1. You should try to define a cut-off radius, which inside this radius the wavefuntion takes on the above value. Also, be careful of your radius step-size: use a small enough step size so your solution doesn't blow up.

The 1/r potential term is best handled by recognizing that the value from the analytical solution at r = 0 for the hydrogen atom is 1. You should try to define a cut-off radius, which inside this radius the wavefuntion takes on the above value. Also, be careful of your radius step-size: use a small enough step size so your solution doesn't blow up.

Thanks a lot!

but I do not mean the 1/r potential

i mean the lapalace operator in polar coordinates, 2D

you will find a term 1/r(d/dr)

This term diverge near the origin.

I do not know how to handle it. Moreover, i am interested in the time-dependent S equation, not the time-independent S equation.