Comprehensive Study Notes

The Numerov Method

We solved the Schrödinger equation exactly for the particle in a box and the harmonic oscillator. For many potential-energy functions $V(x)$, the one-particle, one-dimensional Schrödinger equation cannot be solved exactly. This section presents a numerical method (the Numerov method) for computer solution of the one-particle, one-dimensional Schrödinger equation that allows one to get accurate bound-state eigenvalues and eigenfunctions for an arbitrary $V(x)$.

To solve the Schrödinger equation numerically, we deal with a portion of the $x$ axis that includes the classically allowed region and that extends somewhat into the classically forbidden region at each end of the classically allowed region. We divide this portion of the $x$ axis into small intervals, each of length $s$ (Fig. 4.7). The points $x{0}$ and $x{\text {max }}$ are the endpoints of this portion, and $x{n}$ is the endpoint of the $n$th interval. Let $\psi{n-1}, \psi{n}$, and $\psi{n+1}$ denote the values of $\psi$ at the points $x{n}-s, x{n}$, and $x_{n}+s$, respectively (these are the endpoints of adjacent intervals)

\(
\begin{equation}
\psi{n-1} \equiv \psi\left(x{n}-s\right), \quad \psi{n} \equiv \psi\left(x{n}\right), \quad \psi{n+1} \equiv \psi\left(x{n}+s\right) \tag{4.65}
\end{equation}
\)

Don't be confused by the notation. The subscripts $n-1, n$, and $n+1$ do not label different states but rather indicate values of one particular wave function $\psi$ at points on the $x$ axis separated by the interval $s$. The $n$ subscript means $\psi$ is evaluated at the point $x_{n}$. We write the Schrödinger equation $-\left(\hbar^{2} / 2 m\right) \psi^{\prime \prime}+V \psi=E \psi$ as

\(
\begin{equation}
\psi^{\prime \prime}=G \psi, \quad \text { where } G \equiv m \hbar^{-2}[2 V(x)-2 E] \tag{4.66}
\end{equation}
\)

By expanding $\psi\left(x{n}+s\right)$ and $\psi\left(x{n}-s\right)$ in Taylor series involving powers of $s$, adding these two expansions to eliminate odd powers of $s$, using the Schrödinger equation to express $\psi^{\prime \prime}$ and $\psi^{(\mathrm{iv})}$ in terms of $\psi$, and neglecting terms in $s^{6}$ and higher powers of $s$ (an approximation that will be accurate if $s$ is small), one finds that (Prob. 4.43)

\(
\begin{equation}
\psi{n+1} \approx \frac{2 \psi{n}-\psi{n-1}+5 G{n} \psi{n} s^{2} / 6+G{n-1} \psi{n-1} s^{2} / 12}{1-G{n+1} s^{2} / 12} \tag{4.67}
\end{equation}
\)

where $G{n} \equiv G\left(x{n}\right) \equiv m \hbar^{-2}\left[2 V\left(x{n}\right)-2 E\right]$ [Eqs. (4.66) and (4.65)]. Equation (4.67) allows us to calculate $\psi{n+1}$, the value of $\psi$ at point $x{n}+s$, if we know $\psi{n}$ and $\psi{n-1}$, the values of $\psi$ at the preceding two points $x{n}$ and $x_{n}-s$.

How do we use (4.67) to solve the Schrödinger equation? We first guess a value $E{\text {guess }}$ for an energy eigenvalue. We start at a point $x{0}$ well into the left-hand classically forbidden region (Fig. 4.7), where $\psi$ will be very small, and we approximate $\psi$ as zero at this point: $\psi{0} \equiv \psi\left(x{0}\right)=0$. Also, we pick a point $x{\max }$ well into the right-hand classically forbidden region, where $\psi$ will be very small and we shall demand that $\psi\left(x{\max }\right)=0$. We pick a small value for the interval $s$ between successive points, and we take $\psi$ at $x{0}+s$ as some small number, say, $0.0001: \psi{1} \equiv \psi\left(x{1}\right) \equiv \psi\left(x{0}+s\right)=0.0001$. The value of $\psi{1}$ will not make any difference in the eigenvalues found. If 0.001 were used instead of 0.0001 for $\psi{1}$, Eq. (4.67) shows that this would simply multiply all values of $\psi$ at subsequent points by 10 (Prob. 4.41). This would not affect the eigenvalues [see the example after Eq. (3.14)]. The wave function can be normalized after each eigenvalue is found.

Having chosen values for $\psi{0}$ and $\psi{1}$, we then use (4.67) with $n=1$ to calculate $\psi{2} \equiv \psi\left(x{2}\right) \equiv \psi\left(x{1}+s\right)$, where the $G$ values are calculated using $E{\text {guess }}$. Next, (4.67) with $n=2$ gives $\psi{3}$; and so on. We continue until we reach $x{\max }$. If $E{\text {guess }}$ is not equal to or very close to an eigenvalue, $\psi$ will not be quadratically integrable and $\left|\psi\left(x{\max }\right)\right|$ will be very large. If $\psi\left(x{\max }\right)$ is not found to be close to zero, we start again at $x{0}$ and repeat the process using a new $E{\text {guess }}$. The process is repeated until we find an $E{\text {guess }}$ that makes $\psi\left(x{\max }\right)$ very close to zero. $E{\text {guess }}$ is then essentially equal to an eigenvalue. The systematic way to locate the eigenvalues is to count the nodes in the $\psi$ produced by $E{\text {guess }}$. Recall (Section 4.2) that in a one-dimensional problem, the number of interior nodes is 0 for the ground state, 1 for the first excited state, and so on. Let $E{1}, E{2}, E{3}, \ldots$ denote the energies of the ground state, the first excited state, the second excited state, and so on, respectively. If $\psi{\text {guess }}$ contains no nodes between $x{0}$ and $x{\text {max }}$, then $E{\text {guess }}$ is less than or equal to $E_{1}$; if

FIGURE 4.7 $V$ versus $x$ for a one-particle, onedimensional system.

FIGURE 4.8 The number of nodes in a Numerov-method solution as a function of the energy $E_{\text {guess }}$.

$\psi{\text {guess }}$ contains one interior node, then $E{\text {guess }}$ is between $E{1}$ and $E{2}$ (Fig. 4.8). Examples are given later.

Dimensionless Variables

The Numerov method requires that we guess values of $E$. What should be the order of magnitude of our guesses: $10^{-20} \mathrm{~J}, 10^{-15} \mathrm{~J}, \ldots$ ? To answer this question, we reformulate the Schrödinger equation using dimensionless variables, taking the harmonic oscillator as the example.

The harmonic oscillator has $V=\frac{1}{2} k x^{2}$, and the harmonic-oscillator Schrödinger equation contains the three constants $k, m$, and $\hbar$. We seek to find a dimensionless reduced energy $E{r}$ and a dimensionless reduced $x$ coordinate $x{r}$ that are defined by

\(
\begin{equation}
E{r} \equiv E / A, \quad x{r} \equiv x / B \tag{4.68}
\end{equation}
\)

where the constant $A$ is a combination of $k, m$, and $\hbar$ that has dimensions of energy, and $B$ is a combination with dimensions of length. The dimensions of energy are mass $\times$ length $^{2} \times$ time $^{-2}$, which we write as

\(
\begin{equation}
[E]=\mathrm{ML}^{2} \mathrm{~T}^{-2} \tag{4.69}
\end{equation}
\)

where the brackets around $E$ denote its dimensions, and M, L, and T stand for the dimensions mass, length, and time, respectively. The equation $V=\frac{1}{2} k x^{2}$ shows that $k$ has dimensions of energy $\times$ length $^{-2}$, and (4.69) gives $[k]=\mathrm{MT}^{-2}$. The dimensions of $\hbar$ are energy $\times$ time. Thus

\(
\begin{equation}
[m]=\mathrm{M}, \quad[k]=\mathrm{MT}^{-2}, \quad[\hbar]=\mathrm{ML}^{2} \mathrm{~T}^{-1} \tag{4.70}
\end{equation}
\)

The dimensions of $A$ and $B$ in (4.68) are energy and length, respectively, so

\(
\begin{equation}
[A]=\mathrm{ML}^{2} \mathrm{~T}^{-2}, \quad[B]=\mathrm{L} \tag{4.71}
\end{equation}
\)

Let $A=m^{a} k^{b} \hbar^{c}$, where $a, b$, and $c$ are powers that are determined by the requirement that the dimensions of $A$ must be $\mathrm{ML}^{2} \mathrm{~T}^{-2}$. We have

\(
\begin{equation}
[A]=\left[m^{a} k^{b} \hbar^{c}\right]=\mathbf{M}^{a}\left(\mathrm{MT}^{-2}\right)^{b}\left(\mathrm{ML}^{2} \mathrm{~T}^{-1}\right)^{c}=\mathbf{M}^{a+b+c} \mathrm{~L}^{2 c} \mathrm{~T}^{-2 b-c} \tag{4.72}
\end{equation}
\)

Equating the exponents of each of $\mathrm{M}, \mathrm{L}$, and T in (4.71) and (4.72), we have

\(
a+b+c=1, \quad 2 c=2, \quad-2 b-c=-2
\)

Solving these equations, we get $c=1, b=\frac{1}{2}, a=-\frac{1}{2}$. Therefore,

\(
\begin{equation}
A=m^{-1 / 2} k^{1 / 2} \hbar \tag{4.73}
\end{equation}
\)

Let $B=m^{d} k^{e} \hbar^{f}$. The same dimensional-analysis procedure that gave (4.73) gives (Prob. 4.44)

\(
\begin{equation}
B=m^{-1 / 4} k^{-1 / 4} \hbar^{1 / 2} \tag{4.74}
\end{equation}
\)

From (4.68), (4.73), and (4.74), the reduced variables for the harmonic oscillator are

\(
\begin{equation}
E{r}=E / m^{-1 / 2} k^{1 / 2} \hbar, \quad x{r}=x / m^{-1 / 4} k^{-1 / 4} \hbar^{1 / 2} \tag{4.75}
\end{equation}
\)

Using $k^{1 / 2}=2 \pi \nu m^{1 / 2}$ [Eq. (4.23)] to eliminate $k$ from (4.75) and recalling the definition $\alpha \equiv 2 \pi \nu m / \hbar$ [Eq. (4.31)], we have the alternative expressions

\(
\begin{equation}
E{r}=E / h \nu, \quad x{r}=\alpha^{1 / 2} x \tag{4.76}
\end{equation}
\)

Similar to the equation $E{r} \equiv E / A$ [Eq. (4.68)], we define the reduced potential energy function $V{r}$ as

\(
\begin{equation}
V_{r} \equiv V / A \tag{4.77}
\end{equation}
\)

Since $|\psi(x)|^{2} d x$ is a probability, and probabilities are dimensionless, the normalized $\psi(x)$ must have the dimensions of length ${ }^{-1 / 2}$. We therefore define a reduced normalized wave function $\psi_{r}$ that is dimensionless. From (4.71), $B$ has dimensions of length, so $B^{-1 / 2}$ has units of length ${ }^{-1 / 2}$. Therefore,

\(
\begin{equation}
\psi_{r}=\psi / B^{-1 / 2} \tag{4.78}
\end{equation}
\)

$\psi{r}$ satisfies $\int{-\infty}^{\infty}\left|\psi{r}\right|^{2} d x{r}=1$; this follows from (4.68), (4.78), and $\int{-\infty}^{\infty}|\psi|^{2} d x=1$.
We now rewrite the Schrödinger equation in terms of the reduced variables $x{r}, \psi{r}, V{r}$, and $E_{r}$. We have

\(
\begin{align}
\frac{d^{2} \psi}{d x^{2}}=\frac{d^{2}}{d x^{2}} B^{-1 / 2} \psi{r}=B^{-1 / 2} \frac{d}{d x} \frac{d \psi{r}}{d x} & =B^{-1 / 2} \frac{d}{d x} \frac{d \psi{r}}{d x{r}} \frac{d x{r}}{d x}=B^{-1 / 2} \frac{d\left(d \psi{r} / d x{r}\right)}{d x{r}} \frac{d x{r}}{d x} \frac{d x{r}}{d x} \
\frac{d^{2} \psi}{d x^{2}} & =B^{-5 / 2} \frac{d^{2} \psi{r}}{d x{r}^{2}} \tag{4.79}
\end{align}
\)

since $d x_{r} / d x=B^{-1}$ [Eq. (4.68)]. Substitution of (4.68), (4.77), and (4.79) into the Schrödinger equation $-\left(\hbar^{2} / 2 m\right)\left(d^{2} \psi / d x^{2}\right)+V \psi=E \psi$ gives

\(
\begin{align}
-\frac{\hbar^{2}}{2 m} B^{-5 / 2} \frac{d^{2} \psi{r}}{d x{r}^{2}}+A V{r} B^{-1 / 2} \psi{r} & =A E{r} B^{-1 / 2} \psi{r} \
-\frac{\hbar^{2}}{2 m} \frac{1}{A B^{2}} \frac{d^{2} \psi{r}}{d x{r}^{2}}+V{r} \psi{r} & =E{r} \psi{r} \tag{4.80}
\end{align}
\)

From (4.73) and (4.74), we get $A B^{2}=\hbar^{2} / m$, so $\hbar^{2} / m A B^{2}=1$ for the harmonic oscillator.

More generally, let $V$ contain a single parameter $c$ that is not dimensionless. For example, we might have $V=c x^{4}$ or $V=c x^{2}\left(1+0.05 m^{1 / 2} c^{1 / 2} \hbar^{-1} x^{2}\right)$. (Note that $m^{1 / 2} c^{1 / 2} \hbar^{-1} x^{2}$ is dimensionless, as it must be, since 1 is dimensionless.) The quantity $A B^{2}$ in (4.80) must have the form $A B^{2}=\hbar^{r} m^{s} c^{t}$, where $r, s$, and $t$ are certain powers. Since the term $V{r} \psi{r}$ in (4.80) is dimensionless, the first term is dimensionless. Therefore, $\hbar^{2} / m A B^{2}$ is dimensionless and $A B^{2}$ has the same dimensions as $\hbar^{2} / m$; so $r=2, s=-1$, and $t=0$. With $A B^{2}=\hbar^{2} / m$, Eq. (4.80) gives as the dimensionless Schrödinger equation

\(
\begin{gather}
\frac{d^{2} \psi{r}}{d x{r}^{2}}=\left(2 V{r}-2 E{r}\right) \psi{r} \tag{4.81}\
\psi{r}^{\prime \prime}=G{r} \psi{r}, \quad \text { where } \quad G{r} \equiv 2 V{r}-2 E_{r} \tag{4.82}
\end{gather}
\)

For the harmonic oscillator, $V{r} \equiv V / A=\frac{1}{2} k x^{2} / m^{-1 / 2} k^{1 / 2} \hbar=\frac{1}{2} x{r}^{2}$ [Eqs. (4.73) and (4.75)]:

\(
\begin{equation}
V{r}=\frac{1}{2} x{r}^{2} \tag{4.83}
\end{equation}
\)

Having reduced the harmonic-oscillator Schrödinger equation to the form (4.81) involving only dimensionless quantities, we can expect that the lowest energy eigenvalues will be of the order of magnitude 1 .

The reduced harmonic-oscillator Schrödinger equation (4.82) has the same form as (4.66), so we can use the Numerov formula (4.67) with $\psi, G$, and $s$ replaced by $\psi{r}, G{r}$, and $s{r}$, respectively, where, similar to (4.68), $s{r} \equiv s / B$.

Once numerical values of the reduced energy $E_{r}$ have been found, the energies $E$ are found from (4.75) or (4.76).

Choice of $x{r, 0}, x{r, \text { max }}$, and $s_{r}$

We now need to choose initial and final values of $x{r}$ and the value of the interval $s{r}$ between adjacent points. Suppose we want to find all the harmonic-oscillator eigenvalues and eigenfunctions with $E{r} \leq 5$. We start the solution in the left-hand classically forbidden region, so we first locate the classically forbidden regions for $E{r}=5$. The boundaries between the classically allowed and forbidden regions are where $E{r}=V{r}$. From (4.83), $V{r}=\frac{1}{2} x{r}^{2}$. Thus $E{r}=V{r}$ becomes $5=\frac{1}{2} x{r}^{2}$ and the classically allowed region for $E{r}=5$ is from $x{r}=-(10)^{1 / 2}=-3.16$ to +3.16 . For $E{r}<5$, the classically allowed region is smaller. We want to start the solution at a point well into the left-hand classically forbidden region, where $\psi$ is very small, and we want to end the solution well into the right-hand classically forbidden region. The left-hand classically forbidden region ends at $x{r}=-3.16$ for $E{r}=5$, and a reasonable choice is to start at $x{r}=-5$. [Starting too far into the classically forbidden region can sometimes lead to trouble (see the following), so some trial-and-error might be needed in picking the starting point.] Since $V$ is symmetrical, we shall end the solution at $x{r}=5$.

For reasonable accuracy, one usually needs a minimum of 100 points, so we shall take $s{r}=0.1$ to give us 100 points. As is evident from the derivation of the Numerov method, $s{r}$ must be small. A reasonable rule might be to have $s_{r}$ no greater than 0.1.

If, as is often true, $V \rightarrow \infty$ as $x \rightarrow \pm \infty$, then starting too far into the classically forbidden region can make the denominator $1-G{n+1} s^{2} / 12$ in the Numerov formula (4.67) negative. We have $G{r}=2 V{r}-2 E{r}$, and if we start at a point $x{0}$ where $V{r}$ is extremely large, $G{r}$ at that point might be large enough to make the Numerov denominator negative. The method will then fail to work. We are taking $\psi{0}$ as zero and $\psi{1}$ as a positive number. The Numerov formula (4.67) shows that if the denominator is negative, then $\psi{2}$ will be negative, and we will have produce a spurious node in $\psi$ between $x{1}$ and $x{2}$. To avoid this problem, we can decrease either the step size $s{r}$ or $x{r, \text { max }}-x_{r, 0}$ (see Prob. 4.46).

Computer Program for the Numerov Method

Table 4.1 contains a C ++ computer program that applies the Numerov method to the harmonic-oscillator Schrödinger equation. The fifth and sixth lines declare variables as either integers or double precision. m is the number of intervals between $x{r, 0}$ and $x{r, \text { max }}$ and equals $\left(x{r, \text { max }}-x{r, 0}\right) / s$. cout and cin provide for output to the computer screen and input to the variables of the program. Note the colon at the end of line 13 and the semicolons after most other lines. The three lines beginning with $g[0]=, g[1]=$, and $g[i+1]=$ contain two times the potential-energy function. These lines must be modified if the problem is not the harmonic oscillator. If there is a node between two successive values of $x{r}$, then the $\psi{r}$ values at these two points will have opposite signs (see Prob. 4.45) and the statement $n n=n n+1$ will increase the nodes counter $n n$ by 1 .

You can download a free $\mathrm{C}++$ compiler and integrated development environmentsee the Wikipedia article, List of Compilers. Even simpler, you can enter and run C ++ programs online without downloading anything. The website ideone.com provides this service (after you register) for $\mathrm{C}++$ and other languages. However, at this site, you must enter your complete input into the input area before running the program, since the website does not accept input once the program begins; each input number is separated from its neighbors by a space.

TABLE 4.1 C++ Program for Numerov Solution of the One-Dimensional Schrödinger Equation

#include <iostream>

#include <math.h>

using namespace std;

int main() {<br />
    int m, i, nn;<br />
    double x<span class="hljs-string">[1000]</span>, g<span class="hljs-string">[1000]</span>, p<span class="hljs-string">[1000]</span>, E, s, ss;<br />
    cout << <span class="hljs-string">"Enter initial xr "</span>;<br />
    cin>>x<span class="hljs-string">[0]</span>;<br />
    cout << <span class="hljs-string">"Enter the increment sr "</span>;<br />
    cin>>S;<br />
    cout << <span class="hljs-string">"Enter the number of intervals m "</span>;<br />
    cin>>m;<br />
    label1:<br />
    cout << <span class="hljs-string">"Enter the reduced energy Er (enter 1e10 to quit) "</span>;<br />
    cin>>E;<br />
    if (E > 1e9) {<br />
        cout << <span class="hljs-string">"Quitting"</span>;<br />
        return <span class="hljs-number">0</span>;<br />
    }

    nn=0; p[0]=0;

    p[1]=0.0001;

    x[1]=x[0]+5;

    g[0]=x[0]*x[0]-2*E;

    g[1]=x[1]*x[1]-2*E;

    ss=s*s/12;

    for (i=1; i<=m-1; i=i+1) {<br />
        x<span class="hljs-string">[i+1]</span>=x<span class="hljs-string">[i]</span>+s;<br />
        g<span class="hljs-string">[i+1]</span>=x<span class="hljs-string">[i+1]</span>*x<span class="hljs-string">[i+1]</span>-<span class="hljs-number">2</span>*E;<br />
        p<span class="hljs-string">[i+1]</span>=(-p<span class="hljs-string">[i-1]</span>+<span class="hljs-number">2</span>*p<span class="hljs-string">[i]</span>+<span class="hljs-number">10</span>*g<span class="hljs-string">[i]</span>*p<span class="hljs-string">[i]</span>*ss+g<span class="hljs-string">[i-1]</span>*p<span class="hljs-string">[i-1]</span>*ss)/(<span class="hljs-number">1</span>-g<span class="hljs-string">[i+1]</span>*ss);<br />
        if(p<span class="hljs-string">[i+1]</span>*p<span class="hljs-string">[i]</span><<span class="hljs-number">0</span>)<br />
            nn=nn+<span class="hljs-number">1</span>;<br />
    }

    cout << " Er = " << E << " Nodes = " << nn << " Psir(xm) = " << p[m] << endl;

    goto label1;

}

 

For example, suppose we want the harmonic-oscillator ground-state energy. The program of Table 4.1, with $s{r}=0.1, x{r, 0}=-5$, and $\mathrm{m}=100$ gives the following results. The guess $E{r}=0$ gives a wave function with zero nodes, $\mathrm{nn}=0$, telling us (Fig. 4.8) that the ground-state energy $E{r, 1}$ is above 0 . (Also, the wave function at the rightmost point is found to be $9.94 \times 10^{6}$, very far from 0 . If we now guess 0.9 for $E{r}$, we get a function with one
node, so (Fig. 4.8) 0.9 is between $E{r, 1}$ and $E{r, 2}$. Hence the ground state $E{r}$ is between 0 and 0.9 . Averaging these, we try 0.45 . This value gives a function with no nodes, and so 0.45 is below $E{r, 1}$. Averaging 0.45 and 0.9 , we get 0.675 , which is found to give one node and so is too high. Averaging 0.675 and 0.45 , we try 0.5625 , which gives one node and is too high. We next try 0.50625 , and so on. The program's results show that as we get closer to the true $E{r, 1}, \psi_{r}(5)$ comes closer to zero.

Use of a Spreadsheet to Solve the One-Dimensional Schrödinger Equation

An alternative to a Numerov-method computer program is a spreadsheet.
The following directions for the Excel 2010 spreadsheet apply the Numerov method to solve the harmonic-oscillator Schrödinger equation. (Other versions of Excel can be used with modified directions.)

The columns in the spreadsheet are labeled $\mathrm{A}, \mathrm{B}, \mathrm{C}, \ldots$ and the rows are labeled 1 , $2,3, \ldots$ (see Fig. 4.9 later in this section). A cell's location is specified by its row and column. For example, the cell at the upper left is cell A1. To enter something into a cell, you first select that cell, either by moving the mouse pointer over the desired cell and then clicking the (left) mouse button, or by using the arrow keys to move from the currently selected cell (which has a heavy outline) to the desired cell. After a cell has been selected, type the entry for that cell and press Enter or one of the four arrow keys.

To begin, enter a title in cell A1. Then enter $\mathrm{Er}=$ in cell A3. We shall enter our guesses for $E{r}$ in cell B3. We shall look first for the ground-state (lowest) eigenvalue, pretending that we don't know the answer. The minimum value of $V(x)$ for the harmonic oscillator is zero, so $E{r}$ cannot be negative. We shall take zero as our initial guess for $E_{r}$, so enter 0 in cell B3. Enter $\mathrm{Sr}=$ in cell C3.

Enter 0.1 (the $s{r}$ value chosen earlier) in cell D3. Enter xr in cell A5, Gr in cell B5, and psir in C5. (These entries are labels for the data columns that we shall construct.) Enter -5 (the starting value for $x{r}$ ) in cell A7. Enter $=\mathrm{A} 7+\\( D $\$ 3$ in cell A8. The equal sign at the beginning of the entry tells the spreadsheet that a formula is being entered. This formula tells the spreadsheet to add the numbers in cells A7 and D3. The reason for the $\\) signs in D3 will be explained shortly. When you type a formula, you will see it displayed in the formula bar above the spreadsheet. When you press Enter, the value -4.9 is displayed in cell A8. This is the sum of cells A7 and D3. (If you see a different value in A8 or get an error message, you probably mistyped the formula. To correct the formula, select cell A8 and click in the formula displayed in the formula bar to make the correction.)

Select cell A8 and then click the Home tab and in the Clipboard group click the Copy icon. This copies the contents of cell A8 to a storage area called the Clipboard. Then in the Editing group click the Find \& Select icon and choose Go To. In the Go To box, enter A9:A107 under Reference: and click OK. This selects cells A9 through A107. Then in the Clipboard group click the Paste icon. This pastes the cell A8 formula (which was stored on the Clipboard) into each of cells A9 through A107. To see how this works, click on cell A9. You will see the formula $=A 8+\$ D \$ 3$ in the formula bar. Note that when the cell A8 formula $=A 7+\$ D \$ 3$ was copied one cell down to cell A9, the A7 in the formula was changed to A8. However, the $\\( signs prevented cell D3 in the formula from being changed when it was copied. In a spreadsheet formula, a cell address without $\\) signs is called a relative reference, whereas a cell address with $\$$ signs is called an absolute reference. When a relative reference is copied to the next row in a column, the row number is increased by one; when it is copied two rows below the original row, the row number is increased by two; and so on. Click in some of the other cells in column A to see their formulas. The net result of this copy-and-paste procedure is to fill
the column-A cells with numbers from -5 to 5 in increments of 0.1. (Spreadsheets have faster ways to accomplish this than using Copy and Paste.)

We next fill in the $G{r}$ column. From Eqs. (4.82) and (4.83), $G{r}=x{r}^{2}-2 E{r}$ for the harmonic oscillator. We therefore enter $=A 7 \wedge 2-2 \$ B \$ 3$ in cell B7 (which will contain the value of $G{r}$ at $x{r}=-5$ ). Cell A7 contains the $x_{r}=-5$ value, and the $\wedge$ symbol denotes exponentiation. The denotes multiplication. Cell B3 contains the $E{r}$ value. Next, the rest of the $G{r}$ values are calculated. Select cell B7. Then click Copy in the Clipboard group. Now select cells B8 through B107 by using Find \& Select as before. Then click Paste in the Clipboard group. This fills the cells with the appropriate $G_{r}$ values. (Click on cell B8 or B9 and see how its formula compares with that of B7.)

We now go to the $\psi{r}$ values. Enter 0 in cell C7. Cell C7 contains the value of $\psi{r}$ at $x{r}=-5$. Since this point is well into the classically forbidden region, $\psi{r}$ will be very small here, and we can approximate it as zero. Cell C8 contains the value of $\psi{r}$ at $x{r}=-4.9$. This value will be very small, and we can enter any small number in C 8 without affecting the eigenvalues. Enter $1 \mathrm{E}-4$ in cell C 8 (where E denotes the power of 10 ). Now that we have values in cells C7 and C8 for $\psi{r}$ at the first two points -5.0 and -4.9 [points $x{n-1}$ and $x{n}$ in (4.67)], we use (4.67) to calculate $\psi{r}$ at $x{r}=-4.8$ (point $x{n+1}$ ). Therefore, enter the Eq. (4.67) formula

=(2*C8-C7+5*B8*C8*$D$3^2/6+B7*C7*$D$3^2/12)/(1-B9*$D$3^2/12)

 

in cell C9. After the Enter key is pressed, the value 0.000224 for $\psi{r}$ at $x{r}=-4.8$ appears in cell C9. Select cell C9. Click Copy. Select C10 through C107. Click Paste. Cells C10 through C107 will now be filled with their appropriate $\psi_{r}$ values. As a further check that you entered the C9 formula correctly, verify that cell C10 contains 0.000401 .

Since the number of nodes tells us which eigenvalues our energy guess is between (Fig. 4.8), we want to count and display the number of nodes in $\psi{r}$. To do this, enter into cell D 9 the formula $=\operatorname{IF}(\mathrm{C} 9 * \mathrm{C} 8<0,1,0)$. This formula enters the number 1 into D 9 if C 9 times C 8 is negative and enters 0 into D 9 if C 9 times C 8 is not negative. If there is a node between the $x{r}$ values in A8 and A9, then the $\psi{r}$ values in C8 and C9 will have opposite signs, and the value 1 will be entered into D9. Use Copy and Paste to copy the cell D9 formula into cells D10 through D107. Enter nodes = into cell E2. Enter =SUM(D9:D107) into F2. This formula gives the sum of the values in D9 through D107 and thus gives the number of interior nodes in $\psi{r}$.

Next, we graph $\psi{r}$ versus $x{r}$. Select cells A7 through A107 and C7 through C107 by clicking Find \& Select, clicking Go To, typing A7:A107,C7:C107 in the Reference box, and clicking OK. Then click the Insert tab and in the Chart group click the Scatter icon and then click the chart showing smoothed lines with markers (Scatter with Smooth Lines and Markers). On the chart that appears, right-click Series1 and choose Delete. Right-click the horizontal grid line at $4 \times 10^{6}$ and choose Delete. The spreadsheet will look like Fig. 4.9.

Since $x{r}=5$ is well into the right-hand classically forbidden region, $\psi{r}$ should be very close to zero at this point. However, the graph shows that for our choice $E{r}=0$, the wave function $\psi{r}$ is very large at $x{r}=5$. Therefore, $E{r}=0$ does not give a well-behaved $\psi{r}$, and we must try a different $E{r}$. Cell F2 has a zero, so this $\psi{r}$ has no nodes. Therefore (Fig. 4.8), the guess $E{r}=0$ is less than the true ground-state energy. Let us try $E{r}=2$. Select cell B3 and enter 2 into it. After you press Enter, the spreadsheet will recalculate every column B and column C cell whose value depends on $E{r}$ (all column B and C cells except C 7 and C 8 change) and will then replot the graph. The graph for $E{r}=2$ goes to a very large positive $\psi{r}$ value at $x{r}=5$. Also, cell F2 tells us that $\psi{r}$ for $E{r}=2$ contains two nodes. These are not readily visible on the graph, but the column-C data show that $\psi{r}$ changes sign between the $x_{r}$ values -0.4 and -0.3 and between 1.2 and 1.3 . We are

FIGURE 4.9 Spreadsheet for Numerov-method solution of the harmonic oscillator.

looking for the ground-state eigenfunction, which does not contain a node, or rather, in our approximation, will contain nodes at -5 and 5 . The presence of the two interior nodes shows (Fig. 4.8) that the value 2 for $E{r}$ is not only too high for the ground state, but is higher than $E{r}$ for the first excited state (whose wave function contains one interior node). We therefore need to try a lower value of $E_{r}$.

Before doing so, let us change the graph scale so as to make the nodes more visible. Double click on the $y$ axis of the graph. In the Format Axis dialog box that appears, click Axis Options at the left and click Fixed next to Minimum and Maximum. Replace the original numbers in the Minimum and Maximum boxes with -10 and 10, respectively. Then click Close. The graph will be redrawn with -10 and 10 as the minimum and maximum $y$-axis values, making the two nodes easily visible.

We now change $E{r}$ to a smaller value, say 1.2. Enter 1.2 in cell B3. We get a $\psi{r}$ that goes to a large negative value on the right and that has only one node. The presence of one node tells us that we are now below the energy of the second-lowest state and above the ground-state energy (see Fig. 4.8). We have bracketed the ground-state energy to lie between 0 and 1.2. Let us average these two numbers and try 0.6 as the energy. When we enter 0.6 into cell B3, we get a function with one node, so we are still above the groundstate energy. Since the maximum on the graph is now off scale, it's a good idea to change the graph scale and reset the $y$ maximum and minimum values to 25 and -25 .

We have found the lowest eigenvalue to be between 0 and 0.6 . Averaging these values, we enter 0.3 into B3. This gives a function that has no nodes, so 0.3 is below the groundstate reduced energy, and $E{r}$ is between 0.3 and 0.6 . Averaging these, we enter 0.45 into B3. We get a function that has no nodes, so we are still below the correct eigenvalue. However, if we rescale the $y$ axis suitably (taking, for example, -15 and 30 as the minimum and maximum values), we see a function that for values of $x{r}$ less than 0.2 resembles closely what we expect for the ground state, so we are getting warm. The eigenvalue is now known to be between 0.45 and 0.60 . Averaging these, we enter 0.525 into B3. We get a function with one node, so we are too high. Averaging 0.45 and 0.525 , we try 0.4875 , which we find to be too low.

Continuing in this manner, we get the successive values 0.50625 (too high), 0.496875 (too low), . . , 0.4999995943548 (low), 0.4999996301176 (high). Thus we have found
0.4999996 as the lowest eigenvalue. Since $E_{r}=E / h \nu$, we should have gotten 0.5 . The slight error is due to the approximations of the Numerov method.

Suppose we want the second-lowest eigenvalue. We previously found this eigenvalue to be below 2.0 , so it lies between 0.5 and 2.0. Averaging these numbers, we enter 1.25 into B3. We get a function that has the desired one node but that goes to a large negative value at the right, rather than returning to the $x$ axis at the right. Therefore, 1.25 is too low (Fig. 4.8). Averaging 1.25 and 2.0 , we next try 1.625 . This gives a function with two nodes, so we are too high. Continuing, we get the successive values 1.4375 (low), 1.53125 (high), 1.484375 (low), 1.5078125 (high), and so on.

To test the accuracy of the eigenvalues found, we can repeat the calculation with $s{r}$ half as large and see if the new eigenvalues differ significantly from those found with the larger $s{r}$. Also, we can start further into the classically forbidden region.

Finding eigenvalues as we have done by trial and error is instructive and fun the first few times, but if you have a lot of eigenvalues to find, you can use a faster method. Most spreadsheets have a built-in program that can adjust the value in one cell so as to yield a desired value in a second cell. This tool is called the Solver in Excel. Click the Data tab. In the Analysis group (at the right), see if there is an icon for Solver. (If the Solver icon is missing, click the File tab and click Options at the left. In the Excel Options box, click Add-Ins at the left. In the Manage box select Excel Add-Ins at the bottom and click the Go button. In the Add-Ins box check Solver Add-In and click OK. After a while, the Solver icon will appear in the Analysis group.) To see how the Solver works, enter 0 into cell B3. Click the Solver icon. The Solver Parameters box opens. The $\psi{r}$ value at $x{r}=5$ is in cell C107, and we want to make this value zero. Therefore, in the Set Objective box of the Solver, enter $\$ C \$ 107$. Next to To: click on Value of and enter 0 . We want to adjust the energy so as to satisfy the boundary condition at $x_{r}=5$, so click in the By Changing Variable Cells box and enter $\$ B \$ 3$. Select the Solver Method as GRG Nonlinear. Then click on Solve. When Solver declares it has found a solution, select Accept Solver solution to close the Solver box. The solution found by Excel has 0.500002 in cell B3 (the formula bar will show the full value if you select cell B3). Cell C 107 will have the value -6.38 . Since this value is not close to the desired value of 0 , click the Solver icon and then click Solve, thereby re-running Solver starting from the 0.500002 value. This time Solver gives 0.4999996089 , similar to the value found by hand, and C 107 has the value $-3.66 \times 10^{-9}$. (The Solver uses either the quasi-Newton or the conjugate-gradient method, both of which are discussed in Section 15.10.)

To find higher eigenvalues using the Solver, start with a value in B3 that is well above the previously found eigenvalue. If the program converges to the previous eigenvalue, start with a still-higher value in B3. You can check which eigenvalue you have found by counting the nodes in the wave function. If the Solver fails to find the desired eigenvalue, use trial and error to find an approximate value and then use the Solver starting from the approximate eigenvalue.

The wave function we have found is unnormalized. The normalization constant is given by Eq. (3.93) as $N=\left[\int{-\infty}^{\infty}\left|\psi{r}\right|^{2} d x{r}\right]^{-1 / 2}$. We have $\int{-\infty}^{\infty}\left|\psi{r}\right|^{2} d x{r} \approx \int{-5}^{5} \psi{r}^{2} d x{r} \approx \sum{i=1}^{100} \psi{r, i}^{2} s{r}$, where the $\psi{r, i}$ values are in column B. Enter npsir in E5. In cell D109, enter the formula $=$ SUMSQ(C8:C107)*\$D\$3. The SUMSQ function adds the squares of a series of numbers. Enter $=C 7 / \$ D \$ 109 \wedge 0.5$ in E7. Copy and paste E7 into E8 through E107. Column E will contain the normalized $\psi{r}$ values if $E_{\text {guess }}$ is equal to an eigenvalue.

Excel is widely used to do statistical and scientific calculations, but studies of Excel 2007 and its predecessors found many significant errors [B. D. McCullough and D. A. Heiser, Comput. Statist. Data Anal., 52, 4570 (2008)]. For example, McCullough and Heiser state: "It has been noted that Excel Solver has a marked tendency to stop at a point that is not a solution and declare that it has found a solution." Therefore one should
always verify the correctness of the Solver's solution. A study of Excel 2010 found that many of the errors present in earlier versions of Excel were fixed, but some problems remain [G. Mélard, "On the Accuracy of Statistical Procedures in Excel 2010," available at homepages.ulb.ac.be/~gmelard/rech/gmelard_csda24.pdf]. Mélard noted that "The conclusion is that Microsoft did not make an attempt to fix all the errors in Excel, and this point needs to be made strongly." Mélard found that the Solver yielded results with zero significant-figure accuracy in a substantial fraction of test cases.

Use of Mathcad to Solve the One-Dimensional Schrödinger Equation

Several programs classified as computer algebra systems do a wide variety of mathematical procedures, including symbolic integration and differentiation, numerical integration, algebraic manipulations, solving systems of equations, graphing, and matrix computations. Examples of such computer algebra systems are Maple, Mathematica, MATLAB, Mathcad, and LiveMath Maker. The Numerov procedure can be performed using these programs. One nice feature of Mathcad is its ability to produce animations ("movies"). With Mathcad one can create a movie showing how $\psi{r}$ changes as $E{r}$ goes through an eigenvalue.

Summary of Numerov-Method Steps

Problems 4.30-4.38 apply the Numerov method to several other one-dimensional problems. In solving these problems, you need to (a) find combinations of the constants in the problem that will give a dimensionless reduced energy and length [Eq. (4.68)]; (b) convert the Schrödinger equation to dimensionless form and find what $G{r}\left(x{r}\right)$ in (4.82) is; (c) decide on a maximum $E{r}$ value, below which you will find all eigenvalues; (d) locate the boundaries between the classically allowed and forbidden regions for this maximum $E{r}$ value and choose $x{r, 0}$ and $x{r, \text { max }}$ values in the classically forbidden regions (for the particle in a box with infinitely high walls, use $x{0}=0$ and $x{\max }=l$ ); and (e) decide on a value for the interval $s_{r}$.