Find eigen energies of time-independent Schrödinger equation

You asked for alternative approaches to what you did, so here is one:

A completely different approach to the one-dimensional time-independent Schrödinger equation would be to use matrix techniques. The idea is to eliminate the need for NDSolve entirely. For bound-state problems, you can do this by choosing a basis satisfying the condition of vanishing wave function at infinity and expanding all quantities in that basis.

A particularly nice way of doing this has been pointed out in:

H. J. Korsch and M Glück, "Computing quantum eigenvalues made easy" Eur. J. Phys. 23 (2002) 413-426.

It uses the harmonic-oscillator basis, but also introduces the relation between position and momentum operators on one hand, and the harmonic-oscillator operators on the other hand. This makes it very useful as a teaching method.

I implemented that approach here, with some refinements to satisfy given tolerances and make plots:

---

The potential energy $U(x)$ is assumed to be a polynomial of finite degree in a single variable. Since we're using a harmonic-oscillator basis to find the eigenstates of the Hamiltonian, a finite-degree polynomial translates to a set of finite powers of the position matrix xOp. By increasing the dimension of these matrices while keeping the powers fixed to their finite values, we can achieve convergence of the energy eigenvalues and eigenvectors near the potential minimum.

The larger the matrix dimension, the more of the low-lying eigenvalues will converge to the the true eigenvalues.

First I define position and momentum operators in the basis of harmonic-oscillator eigenfunctions, and a matrix representing the kinetic energy. The dimension of the matrices is dim:

Define basic matrices

xOp[dim_] := 
 SparseArray[{Band[{1, 2}] -> #, Band[{2, 1}] -> #}, {#, #} &@dim] &[
  Table[Sqrt[i], {i, dim - 1}]/Sqrt[2]]

eKin[dim_] := 
 1/2 SparseArray[{Band[{1, 1}] -> Table[i - 1/2, {i, dim}], 
      Band[{1, 3}] -> #, Band[{3, 1}] -> #}, {#, #} &@dim] &[
  Table[-Sqrt[i (i + 1)/2], {i, dim - 2}]/Sqrt[2]]

Spectrum and eigenfunctions

This is the heart of the computation.

The function called spectrum determines the eigenvalues and eigenvectors of a Hamiltonian with the above kinetic energy and a polynomial-type potential energy term. The number of eigenvalues to be returned is given as the first argument, n, and their maximum error is given by tolerance (the latter can be left out if the default value is sufficient):

spectrum[n_, tolerance_: .00001][potential_, var_] /; 
  PolynomialQ[potential, var] := Module[
  {
   x, h, e, v,
   c = CoefficientList[potential, var],
   min = First@Check[
      Minimize[potential, var],
      Print[
       "No absolute potential minimum found. Make sure the leading 
power has a positive coefficient!"];
      Abort[],
      Minimize::natt
      ]
   },
  FixedPoint[
   {#[[1]] + 1, 
     Eigenvalues[
      h = eKin[#[[1]]] + 
        N[Sum[c[[i]] MatrixPower[xOp[#[[1]]], i - 1], {i, 2, 
            Length[c]}] + (c[[1]] - min) IdentityMatrix[#[[
             1]]]], -n]} &, {n + 1, 2 Reverse@Range[n] - 1}, 
   SameTest -> (Max@Abs[#1[[2]] - #2[[2]]] < tolerance &)];
  {e, v} = Reverse /@ Eigensystem[h, -n];
  {e + min, ψe /@ v}
  ]

The second set of arguments above (potential, var) are the potential energy U(x) and the name of the position variable - usually it would be x.

I had to make a change for compatibility with Mathematica version 10 because it throws an error when you try to apply MatrixPower to a numerically singular matrix with power 0. Position matrices can become numerically singular for large dimension, so I have to treat the zeroth power of xOp separately. This wasn't necessary in version 8. Since that term only determines a constant energy offset, it is lumped together with another constant term that I introduce in order to shift the potential minimum to zero (so that the Eigenvalues are sorted in ascending order, starting with the ground state - this shift by min is undone at the end of the calculation).

Output results in position basis

To get the wave function, we implement the sum over basis states to get the wave functions ψe (here, e doesn't have anything to do with parity - the method works independently of parity symmetry). The output format is defined below to yield a shortened representation that doesn't spit out the entire list of coefficients, only listing their number.

Format[ψe[l_List]] := ψe[
  Row[{"[LeftSkeleton]", Length[l], "[RightSkeleton]"}]]
ψe[l_List][x_] := 1/(Pi^(1/4))Exp[-x^2/ 2] l.(1/Sqrt[2^# #!] HermiteH[#, x] &[Range[Length[l]] - 1])

plotSpectrum[data_, potential_, {var_, min_, max_}] :=
 Show[
  Apply[
    Plot[
      Evaluate[#1 + (Last[#1] - First[#1])/(2 + Length[#1])/
          2 Through[#2[var]]], {var, min, max},
      Prolog -> {Gray, Map[Line[{{min, #}, {max, #}}] &, #1]},
      Frame -> True
      ] &,
    data
    ],
  Plot[potential, {var, min, max}]
  ]

Examples

Check that we get the correct numerical eigenvalues for the harmonic oscillator potential $U(x) = frac{1}{2}x^2$ in dimensionless units:

spectrum[10][1/2 x^2, x]

{{0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5}, {ψe[<<15>>], ψe[<<15>>], ψe[<<15>>], ψe[<<15>>], ψe[<<15>>], ψe[<<15>>], ψe[<<15>>], ψe[<<15>>], ψe[<<15>>], ψe[<<15>>]}}

The first list contains the eigenvalues, the second list has the eigenfunctions in shortened form, ready to be plotted. Note from the short forms that that the number of coefficients used in ψe is 15 here, even though I specified that I want to find only 10 eigenvalues in spectrum. This is because these 10 eigenvalues will not be accurate to the given tolerance unless the Hamiltonian matrix is diagonalized in a larger Hilbert space of 15 basis functions. Using FixedPoint, the procedure iteratively increased the number of basis function until none of the first 10 eigenvalues changed by more than the given tolerance in the next iteration.

Example 2: quartic oscillator

Now we'll look at the eigenstates of a quartic potential, shown here:

Plot[-3 x^2 + x^4, {x, -2, 2}]

This is a double-well potential, and the low-lying eigenstates will be localized except for a small tunneling coupling. The potential energy is still a polynomial in $x$, so the method should also work here.

ψ1 = spectrum[5][-10 x^2 + x^4, x] // Last // Last

ψe[<<72>>]

Plot[Evaluate[ψ1[x]], {x, -5, 5}]

If we plot all eigenstates together, the amplitudes are compressed, but you can discern the distinction between even and odd states:

With[{v = -10 x^2 + x^4, n = 10},
 plotSpectrum[
  spectrum[n][v, x], v, {x, -4, 4}]
 ]

As you can see in the quartic-oscillator example, the number of basis states (72) is larger for more complicated potentials. But for the harmonic oscillator test case, it's no problem to get very high accuracy for highly excited levels:

First[spectrum[100][1/2 x^2, x]]

{0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5, 20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5, 30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5, 40.5, 41.5, 42.5, 43.5, 44.5, 45.5, 46.5, 47.5, 48.5, 49.5, 50.5, 51.5, 52.5, 53.5, 54.5, 55.5, 56.5, 57.5, 58.5, 59.5, 60.5, 61.5, 62.5, 63.5, 64.5, 65.5, 66.5, 67.5, 68.5, 69.5, 70.5, 71.5, 72.5, 73.5, 74.5, 75.5, 76.5, 77.5, 78.5, 79.5, 80.5, 81.5, 82.5, 83.5, 84.5, 85.5, 86.5, 87.5, 88.5, 89.5, 90.5, 91.5, 92.5, 93.5, 94.5, 95.5, 96.5, 97.5, 98.5, 99.5}