How can I evaluate more accurate energy eigenvalues from Schrodinger equation using shooting method?

Question

I am trying to use the "shooting method" for solving Schrodinger's equation for a reasonably arbitrary potential in 1D. But the eigenvalues so evaluated in the case of potentials that do not have hard boundaries are not much accurate compared to analytical results.

I thought that the problem could be solved, by making the spatial grid finer, but changing the spatial grid does not practically have any effect on the eigenvalues. I am not making the energy-grid finer, because the job of fining down to the correct eigenvalue is tackled by the bisection method from SciPy, and the wavefunction is evaluated by solving the concerned IVP by odeint from SciPy, these functions are accurate enough.

Finally, changing the 2nd boundary to make the wavefunction die out at deeper part of classically forbidden region also did not bring a practical improvement in the eigenvalue (Changes found only in 9th or 10th place of decimal but made wavefunctions of lower energy state divergent at endpoints to make things worse).

I cannot find what to modify to obtain more accurate eigenvalues. Boundary condition or stepsize? Did my implementation go wrong, or is it due to rounding errors or other "Python things"?

Example: Morse potential

import numpy as np
from scipy.integrate import odeint
from scipy.optimize import bisect


def V(x, xe=1.0, lam=6.0):
    """Morse potential definition"""
    return lam**2*(np.exp(-2*(x- xe)) - 2*np.exp(-(x - xe))) 


def func(y, x):
    """
    Utility function for returning RHS of the differential equation.
    """
    psi, phi = y # psi=eigenfunction, phi=spatial derivative of psi
    return np.array([phi, -(E - V(x))*psi])


def ivp(f, initial1, initial2, X):
    """Solve an ivp with odeint"""
    y0 = np.array([initial1, initial2])
    return odeint(f, y0, X)[:, 0]


def psiboundval(E1): 
    """
    Find out value of eigenfunction at bound2 for energy E1
    by solving ivp.
    """
    global E;
    E = E1
    S = ivp(func, bval1, E1, X)
    return S[(len(S)) - 1] - bval2


def shoot(Erange): 
    """
    Find out accurate eigenvalues from approximate ones in
    Erange by bisect.
    """
    global E
    Y = np.array([psiboundval(E) for E in Erange])
    eigval = np.array([bisect(psiboundval, Erange[i], Erange[i + 1])
                       for i in np.where(np.diff(np.signbit(Y)))[0]])
    return eigval


#%% Solution
xe, lam = 1.0, 6.0 # parameters for potential
# Bval, Bval2 = wavefunction values at x = bound1, bound2
bound1, bound2, bval1, bval2 = 0, xe + 15, 0, 0 
X = np.linspace(bound1, bound2, 1000) # region of integration
Erange = np.geomspace(-lam**2, -0.0001, 100) # region of Energy level searching
print("Numerical results:", np.round(shoot(Erange), 4))
print("Analytical results:",
      [-(lam - n - 0.5)**2 for n in range(0, int(np.floor(lam - 0.5) + 1))])

Output

Numerical results: [-30.2483 -20.2432 -12.2361  -6.2318  -2.2343  -0.2438]
Analytical results: [-30.25, -20.25, -12.25, -6.25, -2.25, -0.25]

For higher energy states, accuracy is seen to decrease.It is desirable that accuracy is of at least up to 4th decimal place (if not more), for all states.

Answer 1

Your problem was the lower limit of integration. It should have been $-x_e$ instead of 0, since $x_e$ is the equilibrium point for the potential and not the minimum distance.

After correcting that, you get the following

#%% Solution
xe, lam = 1.0, 6.0 # parameters for potential
xmax = 10
# Bval, Bval2 = wavefunction values at x = bound1, bound2
bound1, bound2, bval1, bval2 = -xe, xmax, 0, 0 
X = np.linspace(bound1, bound2, xmax) # region of integration
Erange = np.geomspace(-lam**2, -0.01, 100) # region of Energy level searching
print("Numerical results:", np.round(shoot(Erange)[:6], 6))
print("Analytical results:",
      [-(lam - n - 0.5)**2 for n in range(6)])

with results

Numerical results: [-30.25     -20.25     -12.25      -6.25      -2.25      -0.240849]
Analytical results: [-30.25, -20.25, -12.25, -6.25, -2.25, -0.25]

I also tried the method for a simple harmonic oscillator and returned the expected eigenvalues.

As mentioned before, I don't think that this method is the best one for the task. You should try a domain discretization method, such as the Finite Element Method or Finite Differences or a variational method. The latter is commonly used in quantum chemistry codes with Gaussian bases.

The two main drawbacks that I see for this method is that you should have a (finely enough discretized) range for your eigenvalues, that might not be known beforehand. Also, I don't see how to generalize the method to higher-dimensional problems.

A scenario where I think it might be useful is when doing a perturbation analysis. In that case a set of (approximated) eigenvalues is available and the same goes for the eigenfunctions.