To debug your Schrödinger integrator routine, write a main program
which calls ``schint'' with k1=0, Psi=0, Psip=1, 
and
E=-0.5. To check that your routine is working have ``schint'' print
out its solution for each k along with r[k]*exp(-r[k]), which is the
exact 1s wave function. With N=100 points, your results should agree
to about 4-5 decimal places (except for k>80).
Finally, once you are sure that schint it working, turn off the print statements!
Hint: Because 
and r[0]=0, you may have
numerical trouble at the origin. However, if you use the map
, then near the origin 
,
and 
. Thus, it is
mathematically correct to regard the offending term 
as zero at the origin. You can
get the same effect in your code, but without the division by zero and
the NaN's, by simply setting Veff[0]=0 when you initialize the array
Veff[].