Make Numerov find the decaying solution

Adam Zettel · September 25, 2026 · written with an agent

y^{''} = y, y (x) = e^{- x} : outward, the error grows as e^{2 x}

Numerov's method integrates $y^{''} = y$ outward from $y (0) = 1$ , $y (h) = e^{- h}$ , aiming at the decaying solution $e^{- x}$ . By $x = 9$ it is off by a factor of 30, and by $x = 18$ by $2 \times 1 0^{9}$ . Rewrite decaying so that it follows $e^{- x}$ to Numerov's own accuracy, within $2.7 \times 1 0^{- 5}$ , on both intervals.

Open in the workbench Download .vivOutputs below are from the site build’s run; Run all runs it again in your browser.

Look inside this program

Atlas
The program’s map: every function and which calls which.
There’s no recorded run of this program yet.

Numerov's method is the standard way to integrate $y^{''} = f (x) y$ , which is the form of the radial Schrödinger equation. For $y^{''} = k^{2} y$ on a grid of spacing $h$ it is the three-term recurrence

(1 - \frac{h ^{2} k ^{2}}{12}) y_{i + 1} = 2 (1 + \frac{5 h ^{2} k ^{2}}{12}) y_{i} - (1 - \frac{h ^{2} k ^{2}}{12}) y_{i - 1},

accurate to fourth order in $h$ . The equation has two solutions, $e^{- x}$ and $e^{x}$ , and the recurrence carries both. The two starting values are never exactly on the decaying one (they differ from the recurrence's own decaying solution by its truncation error), so a little of $e^{x}$ is always mixed in. Going outward, that part grows by $e^{2 x}$ relative to the part you want. In the region where a bound state's wavefunction decays, every outward integration has this problem.

The function must return $y_{0}, \dots, y_{m}$ at $x_{i} = ih$ , normalized to $y_{0} = 1$ . The tests compare it with $e^{- x}$ at every point. Their tolerances come from posts/refs/make-numerov-find-the-decaying-solution.py, which runs the working method in numpy: its largest relative error is $8.87 \times 1 0^{- 6}$ for $h = 0.15$ on $[0, 9]$ , $5.6 \times 1 0^{- 7}$ for $h = 0.075$ , and $1.84 \times 1 0^{- 5}$ for $h = 0.15$ on $[0, 18]$ . The tests allow about 1.5 times those. Halving $h$ divides the error by 15.9, the factor 16 of a fourth-order method.

Your task

Write these 2 functions:

decaying(m, h): The decaying solution of y'' = y at x_i = ih, i = 0, \ldots, m, with y_0 = 1. As written it runs Numerov outward from y_0 = 1, y_1 = e^{-h}.
max_rel_error(y: Arr, h): The largest relative difference between y and e^{-x} on the grid of spacing h.

Make these tests pass 3 of 5 fail

y(0) = 1 (passes)
with h = 0.15 on [0, 9] every point is within 1.3e-5 of e^(-x), where numpy's inward run reaches 8.87e-6 assert_lt: 29.6647 >= 1.3e-05 (fails)
with h = 0.075 on [0, 9] every point is within 8e-7, where numpy reaches 5.58e-7 assert_lt: 2.00539 >= 8e-07 (fails)
with h = 0.15 on [0, 18], where outward is off by 2e9, every point is within 2.7e-5 (numpy reaches 1.84e-5) assert_lt: 1.9478e+09 >= 2.7e-05 (fails)
halving h divides the error by 12 to 20, as a fourth-order method's should (numpy: 15.9) (passes)

Printed

err_short: 29.6647

tolerance_margin: [1.46498, 1.43461, 1.47079]

Your task

Solutions 0