This is my first question on here. I'm trying to numerically solve the Schrödinger equation for the Woods-Saxon Potential and find the energy eigenvalues and eigenfunctions but I am confused about how exactly this should be done.

I've solved some initial value problems in the past using iterative methods such as Runge–Kutta. I've read that Numerov's method is the way to solve Schrödinger's equation but Wikipedia also describes it as an iterative method for initial value problems.

How do I use it to solve an eigenvalue problem?

This confuses me for the following reasons:

  • Wouldn't iteratively solving the DE require knowledge of the energy eigenvalues to use as input to the calculation? I don't know the eigenvalues yet; they're precisely what I'm trying to calculate.
  • If I did that, wouldn't I simply get a unique solution, instead of a family of eigenfunctions and eigenvalues?

I've seen some mention of "tridiagonal matrices" being generated somehow, but am not sure what the elements of that matrix would be or how that applies to the problem. Leandro M. mentioned that "the discretization defines a finite dimensional (matrix) eigenvalue problem". This seems like the correct road I should be going down, but I haven't been able to find anything that explicitly explains this process or how the matrix is constructed. If this is the correct procedure, how is such a matrix constructed?

  • 1
    $\begingroup$ Would Computational Science be a better home for this question? $\endgroup$
    – Qmechanic
    Commented Mar 16, 2015 at 1:18
  • 2
    $\begingroup$ Much better, yes. (and I notice that the question has already accumulated several more reopen votes since you edited it, so it was probably going to get reopened even without me, which is just how the system is supposed to work!) BTW sorry I didn't get to look at it earlier; I wound up being much busier than I thought last night. $\endgroup$
    – David Z
    Commented Mar 20, 2015 at 5:52
  • 1
    $\begingroup$ If you can solve initial value problems, you can then solve boundary value problem using shooting method. $\endgroup$
    – Ruslan
    Commented Mar 20, 2015 at 7:03
  • 2
    $\begingroup$ I haven't been able to find anything that explicitly explains this process or how the matrix is constructed. Then you have not bothered to open a book on numerical PDEs. $\endgroup$
    – Kyle Kanos
    Commented Mar 30, 2015 at 2:15
  • 1
    $\begingroup$ Why not use the variational method instead? $\endgroup$ Commented Apr 6, 2015 at 14:48

3 Answers 3


I'm glad I can finally answer this.

Numerov's method as described on Wikipedia is not how you want to proceed here. To give you an idea of how to proceed, let's start with a simplified version of the method. What I'm going to do is to just naively discretize the differential operator, like so:

$$ \frac{d^2}{dx^2}\psi \approx \frac{\psi(x-d)-2\psi(x) +\psi(x+d)}{d^2} \equiv \frac{\psi_{n-1} - 2 \psi_n + \psi_{n+1}}{d^2}$$

The last equation is just a definition -- I'm treating space as if it's a lattice with points $d$ apart and I'm calling $\psi_n$ the value of the wavefunction on the $n$th point.. Now the Schrödinger equation reads in this notation:

$$-\frac{\hbar^2}{2m}\frac{\psi_{n-1} - 2 \psi_n + \psi_{n+1}}{d^2} + V_n \psi_n = E \psi_n$$

But this is a matrix equation! Let me be more explicit:

\begin{align} -\frac{\hbar^2}{2md^2}&\left(\begin{array}{cccccccc} -2 & 1 & & & \\ 1 & -2 & 1 & & \\ &\ddots&\ddots&\ddots& \\ & & 1 & -2 & 1 \\ & & & 1 & -2 \\ \end{array}\right) \left(\begin{array}{c} \psi_1 \\ \psi_2 \\ \vdots \\ \psi_{N-1} \\ \psi_N \end{array}\right) +\\ &\qquad\qquad\quad\left(\begin{array}{cccccccc} V_1 & & & & \\ & V_2 & & & \\ & & \ddots & & \\ & & & V_{N-1} & \\ & & & & V_N\\ \end{array}\right)\left(\begin{array}{c} \psi_1 \\ \psi_2 \\ \vdots \\ \psi_{N-1} \\ \psi_N \end{array}\right) = E \left(\begin{array}{c} \psi_1 \\ \psi_2 \\ \vdots \\ \psi_{N-1} \\ \psi_N \end{array}\right) \end{align}

Of course, I had to pick some integer $N$ that just corresponds to placing the system in a box that's "big enough" in order to have a finite system.

It is clear now that what you have is a matrix eigenvalue problem of the form

$$A\psi = E \psi$$

and you may proceed to diagonalize it in whatever way you choose. Note that we call $A$ a tridiagonal matrix, for obvious reasons. You might want to take care of the boundary conditions first -- you do that by setting $\psi_1 = \psi_n = 0$ before you construct the matrix, which corresponds to setting the first and last columns to zero. This will give you some spurious eigenfunctions with zero eigenvalue that you can just discard. If you have different boundary conditions, you're out of luck -- I don't know of a way that makes it work in general.

Now you just have to redo the above with the full Numerov's method, which will be a bit more complicated, and you're all set.

This seems to be the way you're looking for but of course it's not the only way to do this. Griffiths describes one called "wag the dog" that consists of the following: you guess an initial value for an energy and compute the wavefunction however you want (RK, for instance). Chances are it won't be an eigenvalue so the wavefunction will blow up at infinity. It might go to $+\infty$ for large $x$ or it might go to $-\infty$. You now slowly vary the energy until the behavior at infinity "flips", that is, until the tail "wags". That will allow you to constrain the value of a single energy eigenvalue and give you the form of the wavefunction out to some maximum size that increases as the chosen $E$ approaches the exact eigenvalue.

  • 4
    $\begingroup$ You'd better set $\psi_0=\psi_{N+1}=0$ instead of $\psi_1=\psi_N=0$ (i.e. place point with $r=0$ at $n=0$). This way you specifically exclude the points which are of no actual use, and in fact your matrix already represents such choice. Then there'll be no spurious solutions, and as a bonus, you could also solve Coulomb potential problems (because the singularity is excluded from the points you solve for). Note also that since the problem is a 3D problem, and you solve for the radial part of the wavefunction, the solution you find should be divided by $r$. $\endgroup$
    – Ruslan
    Commented Mar 20, 2015 at 9:15
  • $\begingroup$ The method described by Griffiths is more commonly known as the shooting method. $\endgroup$ Commented Oct 8, 2023 at 13:01

You can search for eigenvalues using the bisection method.

Priliminaries: To get the eigenvalues from Numerov method you will need to know the wavefunction at the boundaries. Generally this would mean that you need to set the potential to infinity at the boundaries hence putting the wavefunction to zero at those points. For your potential, modify it as follows:

V = infinity if -50fm>r>50fm V = 0 if |r|>8.5fm V = Woods potential otherwise

The real deal: Now write a program to calculate the wavefunction in the above potential using any arbitrary energy in the domain -50fm<r<50fm. Since, you will start at -50fm the wavefunction will be zero there. If the energy(E1) you chose is an eigenvalue, then you will get the eigenvalue to be 0 at the other boundary too. Otherwise the wavefunction is invalid and implies that the energy you used is not an eigenvalue. Tweak the energy a bit(E2) and recalculate. If the wavefunction at the boundary is still non-zero but has changed sign then you crossed an eigenvalue. In other words there is an eigenvalue between E1 and E2. Now, all you need to do is to use the bisection method on the interval (E1,E2) and find the eigenvalue to the required precision.

Here is a python code that finds an eigen value in the given interval (E1,E2) if it exists.

import numpy as np

rmin = 0
rmax = 0

def vradial(r):
    global rmin 
    global rmax 
    if r < rmin or r > rmax: #This is important. >= and <= won't work
                             #as they interfere with the bc on psi
        return np.inf
    elif r > rmin and r < rmin + 2:
        return (r - rmin - 2)**2
    elif r > rmax - 2 and r < rmax:
        return (r - rmax + 2)**2
        return 0

def f(l,E,r):
    """Calculate the f(r) in the Schrodinger equation of the form
       D2(Psi(r)) = f(r)Psi(r)"""
    if r == 0:
        return np.inf
    return l*(l+1)/(r**2) + vradial(r) - E

def psitophi(psi, l, E, r, delta):
    if psi == 0: return 0 #To handle the 0*infinity case of boundary
    return psi*(1-(f(l, E, r)*delta**2)/12)

def phitopsi(phi, l, E, r, delta):
    #if phi == 0: return 0
    return phi/(1-(f(l, E, r)*delta**2)/12)

def calcpsiradial(rmin, rmax, N, psibcmin, psibcmax, l, E):
    psiarray = []
    r = np.linspace(rmin, rmax, N)
    delta = r[1]-r[0]
    #Assume psi(rmin) = psibcmin and psi(rmin+delta) = 1 and then
    #calculate phi0 and phi1
    phi0 = psitophi(psiarray[0], l, E, r[0], delta) #r[0]=rmin
    phi1 = psitophi(psiarray[1], l, E, r[1], delta)

    #Now populate the psiarray for each value of r
    for i in range(2, len(r)):
        phi2 = 2*phi1 - phi0 + (delta**2)*f(l, E, r[i-1])*psiarray[i-1]
        psi = phitopsi(phi2, l, E, r[i], delta)
        phi0 = phi1
        phi1 = phi2

    return r, psiarray

def normalize(delta, psi):
    area = 0
    for i in range(1,len(psi)-1):
        area = area + abs(psi[i])*delta
    for i in range(len(psi)):
        psi[i] = psi[i]/area
    return psi

def locateEvalueInBracket(rmin, rmax, N, psibcmin, psibcmax, l, e1, e2,\
        tol): #Any value of Psi smaller than psi is
    while abs(e2-e1) > tol:        
        r, psi = calcpsiradial(rmin, rmax, N, psibcmin, psibcmax, l, e1)
        psi1 = psi.pop() 
        r, psi = calcpsiradial(rmin, rmax, N, psibcmin, psibcmax, l, e2)
        psi2 = psi.pop()

        if psi1*psi2 < 0:
            emid = e1 + (e2 - e1) * 0.5
            r, psi = calcpsiradial(rmin, rmax, N, psibcmin, psibcmax, l, emid)
            psimid = psi.pop()

            if psimid*psi1 < 0:
                e2 = emid
            elif psimid*psi2 < 0:
                e1 = emid
        elif psi1*psi2 > 0:
            print "There are either no eigenvalues or too many of them in"+\
            " the given energy interval. For e1={} psi={} and e2={} psi=".format(e1, psi1, e2)+\
            return e1,e2
    return e1,e2

def main():
    import sys
    import matplotlib.pyplot as plt
    global rmin
    global rmax

    e1 = float(sys.argv[1])
    e2 = float(sys.argv[2])
    l = int(sys.argv[3])
    if l == 0:
        rmin = -1
        rmin = 0
    tol =1e-5
    rmax = 1
    psibcmin = 0
    psibcmax = 0
    N = 100

    e1, e2 = locateEvalueInBracket(rmin, rmax, N, psibcmin, psibcmax,\
            l, e1, e2, tol)
    fig = plt.figure()
    ax = fig.add_subplot(111)
    r, psi = calcpsiradial(rmin, rmax, N, psibcmin, psibcmax, l, e1)
    ax.plot(r,psi, 'o')
    r, psi = calcpsiradial(rmin, rmax, N, psibcmin, psibcmax, l ,e2)
    ax.plot(r, psi, 'g')
    ax.set_title("l = {} and E_blue = {}, E_green={}".format(l,e1,e2))

if __name__=="__main__":
  • $\begingroup$ Please use sensible indents while copying and pasting the above code. I can't figure out how to get it to render nicely here. I'll send you a copy if you email me. $\endgroup$
    – gautam1168
    Commented Apr 6, 2015 at 14:11

(This is an substantially expanded version of a comment that was deleted for being mostly just a link)

These kinds of problems easily can be solved by a new approach that I (Grant Bunker, Professor of Physics, IIT, [email protected]) have developed, that I call the Phase Method (http://gbxafs.iit.edu/phase-method). It is fully described and explained in my new book. One of the many examples given is the Woods-Saxon potential. Working example code in Python and Mathematica is posted at that web site. The PM does not use discretization methods, which tend to give poor accuracy for discontinuous potentials, especially for states near the continuum.

The PM is more closely related to the Shooting Method (with the same level of accuracy) but does not require any guessing of eigenvalues or adjusting initial conditions. It automatically maps out the eigenvalues, and can provide wave functions and scattering phase shifts also. Instead of searching for convergent wavefunctions, as is the usual procedure, instead it finds bound state energy eigenvalues by rapidly counting the number of states below any specified energy (using something called the phase plot) and using that state counting function in an evolutionary-tree-like process to map out all the energy eigenvalues, preferably in parallel. It is particularly useful when there are nearly degenerate eigenvalues or you have no idea where they may lurk. Typical relative accuracy using extended precision arithmetic is 15-20 decimal digits. Once accurate eigenvalues are found, wavefunctions are easily calculated, perhaps surprisingly, without having to twiddle boundary conditions.

A simple serial Python implementation is given on the linked web site (as is full code in Mathematica) and is pasted below. The potential in this example is for linear v-shape (Abs[x]/2) but you can change it to anything you like, adjusting the range appropriately for the potential and upper energy limit that is selected. The x-range should extend well beyond the classical turning points at that energy, i.e. the points at which the total energy equals the potential energy. Full explanations for the unconventional aspects of this method are given in the book. Have fun!

import time
import numpy as np
import matplotlib.pyplot as plt
import scipy
import itertools
from scipy.integrate import odeint
# Reference: Grant Bunker, book 'Phase Method' 2024
# if transitions in the Phase Plot are rounded in form
# (they should look like step-functions)
# then extend the range (decrease xmin, and increase xmax)
# if LSODA complains that too much accuracy is requested
# then reduce range (increase xmin, and decrease xmax)
# suggested improvements to this code:
# (to achieve some parity with Mathematica version)
# parallel mapping of bsearch
# extended precision computations
# improved graphics and output
#define potential energy function here
def U(x):

#this is Schr\"odinger's equation as two coupled first order eqns    
def shrek(x, state, e):
    y, yprime = state
    dyprime = 2*(U(x)-e)*y
    dy = yprime
    return [dy, dyprime]

# this code plots the "Phase Plot" for a given energy
def plot_phaseplot(ee):
    y0 = [.01, .01]
    xgrid = np.arange(xmin, xmax, 0.01)
    result = odeint(shrek, y0, xgrid, (ee,), tfirst=True)
    phaseplot = np.tanh(result[:,0])
    plt.axhline( )

#this code 'memo-izes' previously computed number-of-bound-states 
#values so they aren't unnecessarily recomputed -> 3x speed-up
#variable 'savednbs' is a python dictionary of key-value pairs
savednbs = {}
def nbs(ee):
    if (ee) in savednbs:
    y0 = [.1, .1]
    xgrid = np.arange(xmin,xmax, 0.01)
    result = odeint(shrek, y0, xgrid, (ee,), tfirst=True)
    t = np.diff(np.sign(np.tanh(result[:,0])))
    transitions = t[t!=0]/2
    return([ee, nboundstates])

def nstates(ee):
    return(nbs (ee)[1])

def bsearchnbs(intrvl):
    mid =(low+high)/2
    slow = nstates(low)
    shigh = nstates(high)
    smid = nstates(mid)
    newintervals = [ ]
    if slow < smid: 
        newintervals += [[low, mid]]
    if smid < shigh: 
        newintervals += [[mid, high]]

def refine(intervals):
# this maps out the whole eigenvalue spectrum in an initial range 
# using a tree-like binary search
# intervals=[[0.,20.]] example starting interval nested list
# the mapping over the list of intervals would best be done 
# in parallel as in mathematica version
    print("starting interval=",intervals)
    nit=40; # max number of iterations; ~30-40 iterations is OK
    while it <= nit:
        intervals = list(itertools. chain(*phaseplot));#flattens list
        if it%10==0:
           print("===>iteration=",it,";  nintervals=",\
           len(intervals),";  interval midpoints=")
           ngooddigits=7 #typical for 64 bit arithmetic
           roundit=lambda numb: np.round(numb,ngooddigits)
           print(" ")
# after executing the above code, the following commands 
# can be entered into Jupyter (for example) to actually get results
(xmin,xmax) = (-80,80)  
intervals=[startinginterval]; #tees up starting interval for search
evals=refine(intervals); #does the search
print('******************','\nfinal eigenvalue estimates:')
print(evals)  #gives result


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.