# Eigenvalues & Eigenfunctions of the 1-d Schrödinger Equation

**nsolve.f90** is a Fortran 90 module solving the 1-d Schrödinger equation
(especially the radial equation for 3-d problems) using the Numerov method.

## Introduction

**nsolve.f90** finds the solutions \( (f, \varepsilon) \) of the 1-d Schrödinger
equation,
\[ -\frac{1}{2} \frac{d^2 f}{d x^2} + v(x) f(x) = \varepsilon f(x)\,, \]
subject to the boundary conditions
\[ \quad f(a) = f_a, \quad f(b) = f_b\,, \]
where \( f_a \) and \( f_b \) are constants.
The function \( f \) is assumed to be defined on the closed interval \( [a,b] \).

For example, the radial equation for a particle in an isotropoic harmonic
potential is
\[ -\frac{\hbar^2}{2m} \frac{d^2 u}{d r^2} + \left[ \frac{\hbar^2 l (l+1)}{2 m r^2}
+ \frac{1}{2} m \omega^2 r^2 \right] u(r) = E u(r) \,,\] where \( u(r) = r
R(r) \) and \( R(r) \) is the radial part of the wavefunction. When lengths
are expressed in units of the oscillator length \( d=\sqrt{\hbar/m\omega} \) and
energies are expressed in units of \(\hbar \omega\), the
equation to be solved becomes
\[ -\frac{1}{2} \frac{d^2 u}{d r^2} + v(r) u(r) = \varepsilon u(r)\,, \] where
\( v(r) = l(l+1)/2 r^2 + r^2/2 \). The domain of this equation is \(
[0,\infty) \), where the point at 0 is determined by the condition \( u(0) = 0 \).

To solve, we truncate the domain to \( [0,b] \), where \( b \) is large
relative to the oscillator length (e.g. \( b = 10 d \)), and then apply the
Numerov method to integrate
the Schrödinger equation, along with the shooting method to find eigenvalues
and eigenfunctions.

## Download

- Main module: nsolve.f90
Example program (harmonic oscillator): ho.f90

Complete .tar.gz file, with automated tests and other examples: nsolve.tar.gz

## Notes

Currently, the code is designed only to search for bound states.

The code uses a high-order numerical derivative to accurately calculate the
matching point for the two-sided Numerov integration. This allows one to
determine energies to many significant digits (8-10).

Another method would be to compute the matrix elements of the Hamiltonian for
these grid points and diagonalize using a Lanczos method. It would be
interesting to this compare to the Numerov integration method.

## To Do

- Use compensated summation to reduce rounding error in the Numerov
integration.

## References

- S. E. Koonin and D. C. Meredith,
*Computational Physics*, Fortran Version, Addison-Wesley (1990).

*Note that the Fortran code linked from the website is not based on the code given in Ref. 1.*