Eigenvalue problem for differential equations in QM I have a very simple question with regard to numerical methods in physics.
I want to solve the eigenvalue problem for a particle moving in an arbitrary potential. Let's take 1D to be concrete. I.e. I want to find $(E,\psi(x))$ satisfying
\begin{align}
\left[-\frac{1}{2}\partial_x^2 + V(x) \right]\psi(x)=E \psi(x).
\end{align}
Now how do I do it exactly? Naively I would implement the following algorithm:
1) Pick some $E$. 
2) I want to find $\psi(x)$ which is normalizable. So I could pick a large $L > 0 $, set $\psi(-L) = \epsilon > 0$ and $\psi'(-L) = \epsilon'>0$ and numerically integrate from there using the Schrodinger equation. 
3) If I encounter a solution which is exponentially small far to the right of the origin, then I say the solution is normalizable (since it is decaying at $|x|\to\infty$), and I accept the pair $(E,\psi(x))$.
4) I increment $E \to E + dE$ and I repeat the process.
In doing so, I should get the spectrum around my starting value of $E$.
Does this algorithm actually work?? It also seems to me like a very uncontrolled way of doing it; I have no idea how accurate the spectrum is going to be. For example, would changing $L, \epsilon, \epsilon'$ make a difference?
The thing is, I know from Sturm-Liouville theory that the spectrum $E$ is going to be discrete (given $V(x)$ satisfying some nice properties). So the spectrum is going to be a set of measure 0 amongst the entire real line that $E$ lives in. This means that I'm almost surely (i.e. with probability 1) never going to get a solution that is normalizable, and whatever solution I try to numerically integrate from my starting point is always going to blow up having integrated far enough to the right.
So, what algorithm do people use to numerically obtain the spectrum and the eigenvalues? How do I also control accuracy of the spectrum generated?
 A: What you are proposing will work, it is essentially what is known as the shooting method for solving the eigenvalue problem. Note that the eigenfunction is defined up to a multiplicative constant so you can just set $\epsilon$=1 and there is only one parameter $\epsilon'$ to vary in order to achieve a properly decaying solution at infinity. The shooting method is easy to program but it has very limited capability.
A much more powerful approach for dealing with this kind of problem is discretizing all functions and operators on a spatial grid {$x_0, x_1, x_2$, ...}. Then the problem boils down to a linear algebra eigenvalue problem for which there are many reliable canned routines.
A: There are a variety of methods for determining the spectrum of a Hamiltonian system in one or more dimensions. For systems with small numbers of degrees of freedom, a direct matrix approach can be taken, often using some variant of the Sinc Discrete Variable Representation of the Hamiltonian. 
Depending on the size of the Sinc-DVR matrix, either direct $O(n^3)$ diagonalization procedures (for small matrices under $8000\times 8000$), or iterative Krylov space or Lanczos methods (for very large matrices) can be used to determine the relevant parts of the eigenspectrum.
Depending on the degree to which quantum-classical correspondence occurs on the relevant portion of phase space, an appropriate change of basis in the DVR method can be used to achieve significant computational savings by only using basis states which are likely to play a role in the eigenstates of the Hamiltonian. This proceeds by taking advantage of the approximate correspondence between the region of classical phase space occupied by a state with energy $\epsilon$ and the region of the von Neumann lattice space which the quantum state occupies, as illustrated in the following video:
http://www.youtube.com/watch?v=zg-uDK5Iekk
As others have mentioned, your method can be made to work, although I'm not sure that it's used particularly often for these types of calculations, since there are more efficient alternatives such as the ones provided above.
