Tutorial 19. Interactive -- Solve 1D Schrödinger equation numerically as a finite difference boundary value problem

Learning objectives

  • Learn the numerical algorithm for solving Schrödinger equation as a 1D boundary value problem
  • Apply the algorithm to solve 1D Harmonic Oscillator


Earlier we learned that Schrödinger equation can be solved numerically using the numerical techniques for partial differential equations (PDEs) or ordinary differential equations (ODEs). We have learned about the Numerov method in our Interactive tutorial , where we solved the 1D Harmonic Oscillator problem as a finite difference initial value problem: we need to start from one grid point, propagate point by point to obtain the wave function. Obtaining the eigenvalues from the Numerov method is also not very convenient.

In this tutorial, we are going to show that 1D Schrödinger equation can also be solved numerically as a finite difference boundary value problem. As you will see, it is more convenient than the Numerov method. We are going to use the algorithm given in the following paper:

Donald Truhlar, “Finite difference boundary value method for solving one-dimensional eigenvalue equations”, Journal of Computational Physics 10, (1972) 123-132

For a 1D Schrödinger equation in atomic unit ($\hbar=1$)

\[-\frac{1}{2\mu}\frac{d^2}{d R^2}\phi_k(R) + V(R)\phi_k(R) = \epsilon_k \phi_k(R)\]

where $\phi_k(R)$ and $\epsilon_k$ are the $k$-th eigenfunction and eigenvalues of the Hamiltonian. The overall idea of this method is still discretizing the 1D space into $N$ grid points and replacing 2nd order derivative (kinetic energy part) with finite difference. The difference from the Numerov method is that, now we find a set of linear equations the wave function values at the grid points need to satisfy: \(\sum_j F_{ij} \phi_{j,k} = \lambda_k^n \phi_{i,k}\)


\[F_{ij} = \frac{1}{\mu}\left[\delta_{i-1,j}(1-\delta_{i,1})\right] + \delta_{ij}\left[-\frac{2}{\mu}+U^h_{ii}\right] + \frac{1}{\mu}\left[\delta_{i+1,j}(1-\delta_{i,N})\right]\]

\(U^h_{ii} = -2h^2 V(R_i)\) \(\lambda_k^h = -2h^2 \epsilon_k\)

Here, $h$ is the discretization size, $i, j$ denote the indices of the grid points.

Then the problem is transformed into a linear algebra problem. We can just build the F matrix and diagonalize it to obtain the eigenvalues and values of the eigenfunctions of the Hamiltonian simultaneously at all grid points. No point-wise propagation anymore.

Try and learn

We will demonstrate how we can solve the 1D harmonic oscillator problem using the finite difference boundary value techniques. For simplicity, we use atomic units. Then $\hbar=1$. Mass unit is electron mass $m_e$=1. The length unit is bohr. The mass of the harmonic oscillator $\mu=0.5$, and angular frequency $\omega=2$.

Run the following code to explore how the grid size can influence the result


  • First, click “Activate” to activate the code block.
  • Once you see the buttons to “run” at the bottom left corner of the code block, click “run” to run the code.
  • Please be patient. Starting the kernel can be slow sometimes.
  • Select the quantum number $v$ to want to visualize the solution to the wave function of that state. The energy of that state will also be printed out in the title of the plot.
  • Drag the slider bar to adjust the grid size (how many discrete points to use in the range [-10,10])

%matplotlib widget
import ipywidgets as widgets
import numpy as np
import matplotlib.pyplot as plt
import time

hbar = 1
class PDE_BV:
    def __init__(self,m,omega, L, npoints=1000):
        self.m = m
        self.omega = omega
        self.L = L
        self.xlower = -L*0.5
        self.xupper = L*0.5
        self.npoints = npoints
        self.x = np.linspace(self.xlower,self.xupper,self.npoints)
        self.h = self.x[1]-self.x[0]
        self.psi = [None]*self.npoints
        self.E = np.zeros(self.npoints)

    def getV(self, x):
        return 0.5*self.m*self.omega**2*x**2

    def getUii(self, i):
        return  -2*(self.h**2)*self.getV(self.x[i])

    def getF(self):
        F = np.zeros([self.npoints,self.npoints])
        for i in range(0,self.npoints):
            F[i,i] = self.getUii(i) - 2/self.m
            if i > 0:
                F[i,i-1] = 1/self.m
            if i < self.npoints-1:
                F[i,i+1] = 1/self.m

        return F

    def diagonalize(self):
        F = self.getF()
        w,v = np.linalg.eig(F)
        return w,v

    def calc(self):
        w,v = self.diagonalize()
        self.E = - w/2.0/self.h**2
        for k in range(0,len(w)):
            self.psi[k] = v[:,k]
            integral = self.h*np.dot(self.psi[k],self.psi[k])
            self.psi[k] = self.psi[k]/integral**0.5

    def plotWFN(self, v):
        plt.plot(self.x,self.psi[v],label=r'$\psi_v(x)$, k = ' + str(v))
        plt.title(r'$v=$'+ str(v) + r', $E_v$=' + '{:10.4f}'.format(self.E[v]))

### Now define the interactive widgets
# Number of discrete points to use in the range
n_options = widgets.IntSlider(

v_options = widgets.IntSlider(
    description=r'state $v$:',

m1 = 0.5
omega1 = 2
L1 = 20

@widgets.interact(n=n_options, v=v_options)
def update(n,v):
    solver = PDE_BV(m1,omega1, L1, npoints=n)

Further reading

This method is used in solving vibrational self-consistent-field (VSCF), which can be used to solve coupled oscillators found in molecular systems. For more about VSCF, please refer to the following papers:

  1. Joel M. Bowman, “Self‐consistent field energies and wavefunctions for coupled oscillators”, J. Chem. Phys. 68, 608-610 (1978) https://doi.org/10.1063/1.435782
  2. Joel M. Bowman, “The self-consistent-field approach to polyatomic vibrations”, Acc. Chem. Res., 19 202-208 (1986) https://doi.org/10.1021/ar00127a002
Written on November 18, 2021