Tutorial 15. Interactive  Solve Schrödinger Eq Numerically with Numerov method
Learning objectives
 Learn the basics about Numerov method
 Try the interactive python code to solve particle in a box problem numerically
 Understand how grid size can influence the accuracy of gridbased methods.
Background
Analytical solution to the Schrödinger equation is only available for a few model systems. Most systems need to be solved with numerical methods. One type of numerical method is the gridbased method.
There is an excellent introduction to gridbased methods in the online opensource textbook LibreTexts. The article by Prof. Andrei Tokmakoff of the University of Chicago summarized the basic concepts about two gridbased methods: the Numerov method and the Discrete Variable Representation (DVR). Please read that online article to learn about the basic equations used in Numerov method.
Intuitively, the Numerov method is based on finite difference approximation, like other gridbased methods for solving PDE and ODE. It gives us an equation that relates the function value at three neighboring grid points. Given the function value at two initial points, we can propagate and get the function values at the rest grid points onebyone. We can propagate from the left boundary or the right boundary and get two numerical solutions (let’s call them $\psi_\mathrm{left}$ and $\psi_\mathrm{right}$), as shown in the videos below.
 Run Numerov method from left:
 Run Numerov method from right:
This interactive tutorial will use the Numerov method to solve the particleinabox problem in 1 dimension. In the next tutorial, we will use Numerov to solve the Harmonic Oscillator.
Solution
In this example, we use atomic units. Then $\hbar=1$. Mass unit is electron mass $m_e$=1. The energy unit is Hartree $E_h$=1. The length unit is bohr.
For simplicity, we assume that we work on a particle in a box problem where the particle mass $m=1$ ($m_e$), box length L=1 (bohr), and the range inside the box is $x\in [0,1]$. Hence, we have
\[\frac{\hbar^2}{2m}\psi''=E\psi\]for $x\in [0,1]$.
Boundary conditions:
\[\psi(0)=0\] \[\psi(1)=0\]For simplicity, assume we already know the energy of the first eigenstate is
\[\frac{n^2h^2}{8mL^2}=\frac{(2\pi\hbar)^2n^2}{8mL^2}=\frac{\pi^2}{2}\](in bohr). Let’s use Numerov to solve the wave function for the first eigenstate. (In our next tutorial for Harmonic Oscillator, we will show how to search for the eigenvalues (energies) if we don’t know them in advance.)
The Schrödinger equation can be simplified as
\[\psi''=\pi^2 \psi\]Based on Numerov method, the Schrödinger equation of the following form
\[\psi''=k^2(x) \psi\]has the finite difference approximation
\[\psi(x_{n+1}) = 2\psi(x_n)\psi(x_{n1})\delta^2k^2(x_n)\psi(x_n)\]and in this case $k^2(x_n)=\pi^2$
Writing the code is straightforward, except for one immediate problem:
Q: Whether we start the propagation from the left or the right boundary, we only know $\psi(x_0)$, but the propagation needs two known points $\psi(x_0)$ and $\psi(x_1)$ to get the 3rd point $\psi(x_2)$. A: We can choose an arbitrary positive value for $\psi(x_1)$, and later we can normalize the wave function.
Try and learn
Run the following code to explore how the grid size can influence the result
Instructions
 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.
 You will see a plot for the numerical solution of the PIB ground state.

Drag the slider bar to adjust the grid size (how many discrete points to use in the range [0,1])
 Then answer the questions below
 Do the numerical solutions obtained from propagating from the left and from the right match with each other? Do they match when the number of grid points is small (say 3)? Do they match when the number of grid points is large (say >50)?
 Have you noticed that the yaxis scale is changing? It’s because the $\psi$ value at the send initial point is guessed, and the numerical solution is not normalized. How can we handle that?
 Read the title of the plot. How does the runtime change with grid size?
Let’s investigate the integral of the probability density and see how to normalize the wave function
To do this, we are going to use numerical integration of the probability density corresponding to the numerically solved wave function. Please look at the lines below “ # Let’s calculate the numerical integral of psi**2”. We use the Numpy function trapz
to do the numerical integration.
Instructions
 Run the following code block in your browser.
 Drag the slider bar to adjust the grid size (how many discrete points to use in the range [0,1]).
 Read the title of the graph and see what the integral of probability density is. Is it 1 (normalized)?
Now, let’s write a better version of Numerov Solver where we normalize the wave functions
Instructions
 Please look for the lines below
# Calculate the integral of probability distribution
and see how the normalization is done.  Drag the slider bar to adjust the grid size (how many discrete points to use in the range [0,1]).
 Read the title of the graph and see what the integral of probability density is. Is it 1 (normalized)?
 The exact analytical solution is also plotted in the figure (black line). Start with the most sparse grid: npoint=3, gradually increase the number of grid points and see how well the numerical result matches the exact solution.
Take home messages
 “Rome wasn’t built in a day.” Even for this simple example of Numerov solution to PIB, we identify problems and solve them one by one to get a good numerical solution.
 For gridbased numerical methods, the grid size matters!
 There is a tradeoff between the accuracy and computational costs for choosing grid size.