On numerically solving the Schrödinger equation I just read a paper 'A pocket calculator determination of energy eigenvalues' by J Killingbeck
(1979).
Link: http://iopscience.iop.org/0305-4470/10/6/001
I have some questions about section 2.


*

*Why $R<0$ implies the $E$ estimate was too high and $R>1$ implies too low? And is it possible that $0<R<1$ remains?

*Why "if the starting point is chosen at $r=h(l+1)$, then $R(hl)$ can be assigned any
finite value without disturbing the calculation at larger $r$ values" ?

*How to guess a trial energy $E$ (especially the ground state energy)?
(Actually I'm interested in a later paper by Killingbeck (1982) 'Finite-difference methods for eigenvalues'(Link: http://iopscience.iop.org/0022-3700/15/6/009), in which the author applies the method in 'A pocket calculator determination of energy eigenvalues' to some potentials with highly singularity at origin)

mod note: copying the text of whoplisp's post here since it isn't actually an answer to the question
They convert the Schroedinger equation into a finite difference equation
(instead of a variational approach):
$-\frac{1}{2}\Delta\psi +V\psi=E\psi$
$V$ is a central potential
Ansatz with solid harmonic $\xi_l$ of degree $l$ (e.g. $x$ or $xy$)
$\psi=\xi_l\phi(r)$
leads to
$2r(V-E)\phi=rD^2\phi+(2l+2)D\phi$
Notation: $D=\frac{d}{dr}$
finite difference approximation
$h^2D^2\phi \rightarrow \phi(r+h)+\phi(r-h)-2\phi(r)$
$2hD\phi \rightarrow \phi(r+h)-\phi(r-h)$
crucial step: introduce ratio variable $R(r)$:
$\phi(r+h)=R(r)\phi(r)$
convert Schroedinger equation in recursive equation to calculate $R(r)$:
$2r[1+h^2(V-E)]=R(r)[r+h(l+1)]+[r-h(l+1)](R(r-h))^{-1}$

the corresponding 1D Schroedinger equation is:
$2[1+h^2(V-E)]=R(x)+(R(x-h))^{-1}$

if $r=h(l+1)$ is chosen then $R(hl)$ can be assigned any finite value without
disturbing the calculation at large $r$ values
to find the ground states with angular momentum $l$ start with an estimate for $E$,
some non-zero $R(hl)$ and some small $h$. $R(r)$ as $r$ increases will either become negative or pass through unity from below.
first thing happens when guessed $E$ was too high
second thing indicates $E$ was too low (because $\psi$ increases with $r$)
you just try various values and 'sandwhich' the solution
example for $V=\lambda r-1/r$:
$\lambda=1$, $h=.1$, state=2p: 1.9759
 A: I can't entirely make sense of some the claims in that letter, but this looks like a version of solving for the bound states of potential which is fairly standard. It's in Griffiths Quantum Mechanics I think. I'll write it down here for you and maybe you can piece together the rest. Anyway it goes like this:
We're trying to find the values of $E$ for which the radial Schrodinger equation 
$-\frac{1}{r^2}\frac{\partial}{\partial r}(r^2\frac{\partial\psi}{\partial r} )+ \frac{\ell(\ell+1)}{r^2}\psi + V(r)\psi = E\psi$
has normalizable solutions. On the other hand that's a second order ODE for any $E$ so it has to have two solution - they're just not going to be normalizable. Specifically they are going to blow up quite badly as $r\rightarrow\infty$. So we are just going to guess a value of $E$, numerically solve that equation by finite differences and see if the solutions blow up at infinity or not. When we found a value of $E$ where they don't blow up we know we found a bound state.
As for the particular issues that come up when we are looking for the ground state. 


*

*$R < 0$ means your wavefunction has a node (it changes sign) and $R > 1$ means your wavefunction is increasing. It is the case that the ground state wavefunction has no nodes so you must be looking in the vicinity of one of the excited states. As for $R > 1$, while asymptotically your wavefunction has to decrease, I don't see why its not possible that it could increase at some point. So this does not seem to be valid termination rule in general (though given a particular potential you might be able to say more).

*See that $R(h(l+1))$ does not depend on $R(hl)$ since the prefactor is zero. After that $R(x)$ depends only on $R(x-h)$ so $R(hl)$ never enters. 

*This method is very computationally cheap so why not just try a lot of energies? You can always bound the energy from above by the standard variational methods. Also just dimensional analysis could give you the order of magnitude.
