Initial Value Problem Methods: A Computational Laboratory

Compare numerical methods for \(y' = f(t,y),\ y(t_0)=y_0\) — Euler, improved Euler, backward Euler, trapezoidal, Adams–Bashforth, Adams–Moulton, BDF, Taylor, and Runge–Kutta. Tick several methods, vary the step size live, and study accuracy, convergence order, and absolute stability on stiff and non-stiff problems.

The Initial Value Problem

A first-order initial value problem (IVP) asks for a function \(y(t)\) satisfying

\[ y'(t) = f(t,y), \qquad y(t_0)=y_0, \qquad t\in[t_0,\,t_{\text{end}}]. \]

A numerical method replaces the continuous solution by a sequence \(y_0, y_1, y_2,\dots\) on a grid \(t_n = t_0 + nh\), where \(h\) is the step size. The recipe for getting \(y_{n+1}\) from earlier values defines the method.

One-step vs. multistep

One-step methods use only \(y_n\) to produce \(y_{n+1}\): Euler, improved Euler (Heun), the Taylor methods, and the Runge–Kutta family. Backward Euler and the trapezoidal rule are implicit one-step methods.
Linear multistep methods use several past values. Adams–Bashforth (AB) are explicit; Adams–Moulton (AM) and the backward differentiation formulas (BDF) are implicit. BDF use only the newest \(f\) value and are the standard methods for stiff problems.
Predictor–corrector methods pair an explicit predictor with an implicit corrector applied once, gaining much of the corrector’s accuracy and stability without an implicit solve: the Adams pairs ABM2/3/4 and the fourth-order Milne–Simpson and Hamming methods.
Explicit methods give \(y_{n+1}\) by a direct formula. Implicit methods define \(y_{n+1}\) through an equation that must be solved each step (here by Newton’s method) — more work per step, but far better stability.

Predictor–corrector methods

A predictor–corrector pair runs an explicit method and an implicit one in tandem, in PECE form: Predict \(y^{*}_{n+1}\) with the explicit formula, Evaluate \(f(t_{n+1},y^{*}_{n+1})\), Correct with the implicit formula using that slope, then Evaluate once more. The corrector is applied a single time rather than iterated to convergence, so the whole step stays explicit — no Newton solve, two \(f\) evaluations per step — yet inherits much of the corrector’s accuracy and its far larger stability region.

The Adams pairs ABM2, ABM3, ABM4 use an Adams–Bashforth predictor with the matching Adams–Moulton corrector; ABM4 is the classic combination.
Milne–Simpson and Hamming are both fourth order and share Milne’s predictor, differing only in the corrector — and that difference is instructive. Milne–Simpson is weakly unstable: a parasitic root of its characteristic polynomial sits on the unit circle, so error oscillations never decay. Hamming’s corrector pushes those parasitic roots strictly inside the circle, curing the instability. The Milne vs Hamming example makes the contrast visible.

Starting a multistep method

A \(k\)-step method needs \(k\) values before its formula can be applied, but the problem supplies only \(y_0\). The missing starting values \(y_1,\dots,y_{k-1}\) are generated here with an explicit Runge–Kutta starter:

Order 2 methods (AB2, ABM2) start with improved Euler (Heun);
Order 3 methods (AB3, AM3, ABM3) start with third-order Runge–Kutta;
Order 4 methods (AB4, AM4, ABM4, Milne–Simpson, Hamming) start with classical RK4.

The starter is chosen to match the multistep method’s order: to keep the global error \(O(h^p)\), the starting values must themselves be accurate to \(O(h^{p})\). A lower-order starter (say improved Euler used for AB4) would pollute those first points and visibly drop the measured convergence rate in the Convergence panel.

The BDF methods are meant for stiff problems, where an explicit starter would itself blow up. Their starting values are therefore generated with Radau IIA (a three-stage implicit Runge–Kutta method of order 5 that is L-stable), so the startup stays stable on stiff problems while its high order keeps the starting error well below the BDF order.

Order of accuracy

A method has order \(p\) if its local truncation error per step is \(O(h^{p+1})\), giving a global error \(O(h^p)\) over a fixed interval. Halving \(h\) then cuts the error by a factor \(2^p\). The Convergence panel measures this slope directly on a log–log plot.

Measuring the error

The global error can be summarized by a single number in several ways; the norm selector controls which is used in the error table and the convergence plot:

\(L^2\) (RMS): \(\sqrt{\frac1N\sum_n e_n^2}\) — the error averaged over the whole trajectory. Robust and gives the cleanest convergence slopes; a good default.
Max (\(L^\infty\)): \(\max_n|e_n|\) — the worst error anywhere; this is the textbook “global error”.
Final-time: \(|y_N-y(t_N)|\) — only the endpoint. Useful when that value is what you want, but it can read accidentally small when the error oscillates (e.g. the flame problem, which sits at \(y=1\) by \(t_{\text{end}}\)).

All consistent norms share the same asymptotic order \(p\); they differ in robustness and in what they emphasize.

Absolute stability

Apply a method to the linear test equation \(y'=\lambda y\) with step \(h\). Each step multiplies the numerical solution by an amplification factor that depends only on \(z=h\lambda\). The method is absolutely stable when that factor has modulus \(\le 1\), so errors do not grow. The set of such \(z\) in the complex plane is the method’s region of absolute stability, drawn in the Stability panel.

Explicit methods (Euler, RK, AB, and the predictor–corrector pairs) have bounded stability regions — on a stiff problem (very negative \(\lambda\)) the step \(h\) must be tiny or the solution blows up. The predictor–corrector regions are larger than their bare predictors but still bounded.
Backward Euler, the trapezoidal rule, and BDF2 are A-stable: their stability regions contain the entire left half-plane, so they stay bounded for any \(h>0\) when \(\operatorname{Re}\lambda<0\). BDF3–5 are A(\(\alpha\))-stable — stable in a wide wedge around the negative real axis. This is why implicit methods are used for stiff problems.

Adaptive step size

Everything above uses a fixed step \(h\). Real solvers instead choose \(h\) automatically to keep the local error near a user-set tolerance, taking small steps through sharp transients and large steps where the solution is smooth. Turn on Adaptive step size to switch the tool into this mode for RK4, Tsitouras 5(4), Dormand–Prince 5(4), Runge–Kutta–Fehlberg 4(5) and DOP853 8(5,3); the step-count control is replaced by a tolerance slider and the other methods are disabled.

Tsit5 and DOPRI5 are embedded pairs: alongside the order-5 step they compute an order-4 estimate, and the difference \(\|y_5-y_4\|\) estimates the local error for free. RK4 has no embedded estimate, so it uses step-doubling — one step of size \(h\) versus two of \(h/2\), with error \(\approx|y_{h/2}-y_h|/15\). DOP853 carries two embedded estimators (5th- and 3rd-order): it blends them as \(\text{err}=|h|\,\hat e_5^2/\sqrt{\hat e_5^2+0.01\,\hat e_3^2}\), the 3rd-order term guarding the 5th against the spurious zeros a lone estimator can hit — and both reuse the twelve stage slopes the order-8 step already evaluated, so the estimate costs no extra work.
After each trial step the error estimate is compared with the tolerance: the step is accepted or rejected, and \(h\) is grown or shrunk by \(0.9\,(\text{tol}/\text{err})^{1/(\hat p+1)}\) — exponent \(1/5\) for the order-4 estimators, \(1/8\) for DOP853's order-7 estimator (growth is held back right after a rejection to keep the controller stable).

Because adaptive stepping has no single \(h\), the tool shows just three panels in this mode. The solution plot marks the non-uniform accepted steps; the error growth panel still tracks the error at those nodes; and the convergence panel becomes a work–precision diagram — achieved error against the number of accepted steps over a sweep of tolerances, so lower-and-to-the-right is more efficient. The stability region panel is hidden, since absolute stability is a fixed-step notion. The results table reports steps accepted, rejected, and the function-evaluation count.

Using this tool

Tick any number of methods in the selector to overlay their approximations on the solution plot.
Pick a built-in example (stiff and non-stiff) — its equation and known exact solution are shown above the solution plot — or enter a custom \(f(t,y)\). Error is measured against the example’s exact solution; a custom problem (or any edited equation) uses a high-accuracy numerical reference instead — dense RK4, automatically switching to an A-stable implicit solver if the problem is stiff. Exact solutions are never typed in, so a wrong formula can’t silently distort the error and convergence panels.
Drag the step-size slider with Live update on to watch accuracy and stability change instantly.
For parameterized examples like \(y'=\lambda y\) (exact \(y_0e^{\lambda t}\)), drag \(\lambda\) from gentle to stiff.
The plot supports scroll-zoom, drag-to-pan, and box-zoom; the graph height is adjustable.

Euler & one-step methods

Forward (Explicit) Eulerexplicitorder 1

The simplest method: follow the tangent line for one step.

\[ y_{n+1} = y_n + h\,f(t_n,y_n) \]

It is the order-1 Taylor method and the first-order Adams–Bashforth method. Local error \(O(h^2)\), global error \(O(h)\). Stability function \(R(z)=1+z\); the stability region is the disc \(|1+z|\le1\), a unit circle centered at \(-1\). Cheap but only conditionally stable, so stiff problems force a tiny \(h\).

Improved Euler / Heun (RK2) explicitorder 2

A predictor–corrector that averages the slope at the start and at the Euler-predicted endpoint — the explicit trapezoidal rule.

\[ k_1=f(t_n,y_n),\quad k_2=f(t_n+h,\;y_n+h k_1),\quad y_{n+1}=y_n+\dfrac{h}{2}(k_1+k_2) \]

Second order, \(R(z)=1+z+\dfrac{z^2}{2}\). Two function evaluations per step buy an extra order over Euler.

Backward (Implicit) Euler implicitorder 1

Evaluate the slope at the new point, giving an equation for \(y_{n+1}\) solved here by Newton’s method.

\[ y_{n+1} = y_n + h\,f(t_{n+1},\,y_{n+1}) \]

First order with \(R(z)=\dfrac{1}{1-z}\). It is A-stable (indeed L-stable): for \(\operatorname{Re}\lambda<0\) it stays bounded for every \(h\), making it a workhorse for stiff problems. It is also the first-order Adams–Moulton method.

Trapezoidal Rule (Crank–Nicolson) implicitorder 2

Average the slopes at both ends — the implicit trapezoidal / second-order Adams–Moulton method.

\[ y_{n+1} = y_n + \dfrac{h}{2}\big[f(t_n,y_n)+f(t_{n+1},y_{n+1})\big] \]

Second order, \(R(z)=\dfrac{1+z/2}{1-z/2}\). A-stable; the boundary is exactly the imaginary axis, so it preserves oscillations well but damps stiff transients only slowly.

Runge–Kutta methods

Midpoint Runge–Kutta (RK2)explicitorder 2

Use the slope at the midpoint of the step.

\[ k_1=f(t_n,y_n),\quad k_2=f\!\big(t_n+\dfrac h2,\,y_n+\dfrac h2 k_1\big),\quad y_{n+1}=y_n+h\,k_2 \]

Second order, same stability function as Heun, \(R(z)=1+z+\dfrac{z^2}{2}\).

Runge–Kutta 3 explicitorder 3

Kutta’s classical three-stage method.

\[ \begin{aligned}k_1&=f(t_n,y_n),\ k_2=f(t_n+\dfrac h2,y_n+\dfrac h2 k_1),\\ k_3&=f(t_n+h,\,y_n-h k_1+2h k_2),\\ y_{n+1}&=y_n+\dfrac h6(k_1+4k_2+k_3).\end{aligned} \]

Third order, \(R(z)=1+z+\dfrac{z^2}{2}+\dfrac{z^3}{6}\).

Runge–Kutta 4 (classical) explicitorder 4

The famous four-stage RK4 — the standard general-purpose explicit method.

\[ \begin{aligned}k_1&=f(t_n,y_n),\quad k_2=f(t_n+\dfrac h2,y_n+\dfrac h2 k_1),\\ k_3&=f(t_n+\dfrac h2,y_n+\dfrac h2 k_2),\quad k_4=f(t_n+h,y_n+h k_3),\\ y_{n+1}&=y_n+\dfrac h6(k_1+2k_2+2k_3+k_4).\end{aligned} \]

Fourth order, \(R(z)=\sum_{j=0}^{4}\dfrac{z^j}{j!}\). Four evaluations per step; excellent accuracy for non-stiff problems.

Dormand–Prince 5(4) explicitorder 5

A seven-stage embedded 5(4) pair — the method behind MATLAB’s ode45 and SciPy’s RK45. The order-5 solution advances the step while the embedded order-4 solution estimates the local error for adaptive step control. It is FSAL (first-same-as-last): the last stage equals the first of the next step, so a step costs six new \(f\) evaluations.

Here the step size is fixed, so only the order-5 formula is used (the error estimate would drive adaptivity). Its stability function is a degree-6 polynomial \(R(z)=1+z+\dots+\dfrac{z^5}{120}+\dfrac{z^6}{600}\) — note the \(z^6\) term is \(\tfrac1{600}\), not \(\tfrac1{720}\); the real stability interval reaches about \(-3.3\).

Tsitouras 5(4) explicitorder 5

Tsitouras’s 2011 seven-stage 5(4) pair — the default Tsit5 integrator in Julia’s OrdinaryDiffEq. Its coefficients are optimized to minimize the order-5 truncation error while keeping a good stability region, and it is also FSAL (six new \(f\) evaluations per fixed step).

Like Dormand–Prince it is used here at fixed step (order-5 formula only). Its degree-6 stability polynomial differs from DOPRI5 in the \(z^6\) coefficient, giving a slightly larger real stability interval of about \(-3.5\).

Runge–Kutta–Fehlberg 4(5) explicitorder 4

Fehlberg’s 1969 embedded pair — the original adaptive Runge–Kutta method. Six \(f\) evaluations yield both a fourth-order solution and a fifth-order one; their difference estimates the local error and drives the step size. In the classic form it advances with the fourth-order formula and uses the fifth only for error control.

\[ \begin{array}{c|cccccc} 0\\ \tfrac14 & \tfrac14\\ \tfrac38 & \tfrac{3}{32} & \tfrac{9}{32}\\ \tfrac{12}{13} & \tfrac{1932}{2197} & -\tfrac{7200}{2197} & \tfrac{7296}{2197}\\ 1 & \tfrac{439}{216} & -8 & \tfrac{3680}{513} & -\tfrac{845}{4104}\\ \tfrac12 & -\tfrac{8}{27} & 2 & -\tfrac{3544}{2565} & \tfrac{1859}{4104} & -\tfrac{11}{40}\\ \hline & \tfrac{25}{216} & 0 & \tfrac{1408}{2565} & \tfrac{2197}{4104} & -\tfrac15 & 0 \end{array} \]

Tick Adaptive step size to let it choose \(h\) automatically to the requested tolerance; at fixed step it is a plain fourth-order method.

Dormand–Prince 8(5,3) (DOP853) explicitorder 8

The famous high-order DOP853 of Hairer, Nørsett & Wanner — a twelve-stage explicit method of order 8, with an embedded 5th- and 3rd-order pair for error control (the “8(5,3)”). It is the workhorse non-stiff integrator behind SciPy’s DOP853 and the original dop853.f.

At fixed step it uses only its order-8 propagation formula. With twelve stages the cost per step is high, but the error drops as \(O(h^8)\): on the convergence panel its slope is the steepest of all the methods, reaching the round-off floor after only a few refinements.

Tick Adaptive step size to run the full DOP853: the embedded 5th- and 3rd-order pair (the “8(5,3)”) now drives error-controlled stepping, exactly as in SciPy’s DOP853 and the reference dop853.f. Because the two estimators reuse the twelve stage slopes, each step still costs just twelve \(f\) evaluations while automatically shrinking \(h\) through sharp transients and stretching it where the solution is smooth.

An order-8 method takes very large steps, so its accepted nodes are sparse. To draw the trajectory faithfully the tool also builds DOP853’s dense output — Hairer & Wanner’s 7th-order continuous interpolant — and samples it between steps. The plotted curve then follows the true solution rather than straight chords between nodes, with the ● markers showing where the accepted steps actually land. The interpolant costs three extra stage evaluations per step and is the same one used for event location and continuous output in production solvers.

Gauss–Legendre (2-stage) implicitorder 4

A fully implicit two-stage Runge–Kutta method whose nodes are the Gauss–Legendre quadrature points \(c=\tfrac12\pm\tfrac{\sqrt3}{6}\). With only two stages it attains order 4 — the maximum possible — and is both A-stable and symplectic.

\[ \begin{array}{c|cc} \tfrac12-\tfrac{\sqrt3}{6} & \tfrac14 & \tfrac14-\tfrac{\sqrt3}{6}\\[2pt] \tfrac12+\tfrac{\sqrt3}{6} & \tfrac14+\tfrac{\sqrt3}{6} & \tfrac14\\[2pt] \hline & \tfrac12 & \tfrac12 \end{array} \]

Because both stages are coupled, each step solves a \(2\times2\) nonlinear system for the stage slopes by Newton iteration. Its stability function is the \((2,2)\) Padé approximant of \(e^z\), \(R(z)=\dfrac{1+z/2+z^2/12}{1-z/2+z^2/12}\); the boundary is exactly the imaginary axis, so the whole left half-plane is stable. Unlike backward Euler it is not L-stable (\(|R|\to1\) as \(z\to-\infty\)), so it preserves oscillations but damps very stiff transients only slowly.

SDIRK 2 implicitorder 2

A two-stage singly diagonally implicit Runge–Kutta method (SDIRK): the diagonal of \(A\) is a single repeated value \(\gamma=1-\tfrac{\sqrt2}{2}\), so the stages are solved one after another (each a scalar implicit equation) instead of as a coupled system. It is order 2 and L-stable.

\[ \begin{array}{c|cc} \gamma & \gamma & 0\\[2pt] 1 & 1-\gamma & \gamma\\[2pt] \hline & 1-\gamma & \gamma \end{array}\qquad \gamma=1-\tfrac{\sqrt2}{2} \]

The bottom row of \(A\) equals the weights \(b\), so the method is stiffly accurate and L-stable (\(R(\infty)=0\)): it damps very stiff transients strongly like backward Euler, but with second-order accuracy — a cheap, robust implicit choice for stiff problems.

Radau IIA (3-stage) implicitorder 5

A fully implicit three-stage Runge–Kutta method whose nodes are the Radau quadrature points (with \(c_3=1\)). With three stages it attains order 5 — the maximum possible — and is L-stable. This is the order-5 core of Hairer & Wanner’s RADAU5, the standard high-order integrator for stiff problems.

\[ \begin{array}{c|ccc} \frac{4-\sqrt6}{10} & \frac{88-7\sqrt6}{360} & \frac{296-169\sqrt6}{1800} & \frac{-2+3\sqrt6}{225}\\[4pt] \frac{4+\sqrt6}{10} & \frac{296+169\sqrt6}{1800} & \frac{88+7\sqrt6}{360} & \frac{-2-3\sqrt6}{225}\\[4pt] 1 & \frac{16-\sqrt6}{36} & \frac{16+\sqrt6}{36} & \frac19\\[4pt] \hline & \frac{16-\sqrt6}{36} & \frac{16+\sqrt6}{36} & \frac19 \end{array} \]

All three stages are coupled, so each step solves a \(3\times3\) nonlinear system for the stage slopes by Newton iteration. The bottom row of \(A\) equals the weights \(b\), so the method is stiffly accurate and L-stable (\(R(\infty)=0\)): unlike Gauss–Legendre it damps very stiff transients strongly, combining high order with strong stiff decay. It is used here to generate the starting values for every BDF method.

Lobatto IIIA (3-stage) implicitorder 4

A fully implicit three-stage Runge–Kutta method whose nodes are the Lobatto quadrature points \(c=0,\tfrac12,1\) — both endpoints are included. With three stages it attains order 4. The two-stage member of the same family is exactly the trapezoidal rule, so Lobatto IIIA is its natural higher-order extension.

\[ \begin{array}{c|ccc} 0 & 0 & 0 & 0\\[2pt] \tfrac12 & \tfrac{5}{24} & \tfrac13 & -\tfrac{1}{24}\\[2pt] 1 & \tfrac16 & \tfrac23 & \tfrac16\\[2pt] \hline & \tfrac16 & \tfrac23 & \tfrac16 \end{array} \]

The first row of \(A\) is zero (so the first stage is just \(f(t_n,y_n)\)) while the remaining stages are coupled, so each step solves a nonlinear system for the stage slopes by Newton iteration. The bottom row of \(A\) equals the weights \(b\), so the method is stiffly accurate. It is A-stable but not L-stable (\(|R(\infty)|=1\)): like Gauss–Legendre it preserves oscillations and damps very stiff transients only slowly.

Taylor series methods

Taylor method, order 2explicitorder 2

Match the Taylor expansion of the solution. The needed total derivatives of \(y\) are obtained by symbolic differentiation of \(f\):

\[ y_{n+1}=y_n+h\,y'_n+\dfrac{h^2}{2}y''_n,\qquad y''=f_t+f_y f. \]

Order 2, same stability function as RK2. Taylor methods need derivatives of \(f\), so they shine when \(f\) is given by a clean formula.

Taylor method, order 3 explicitorder 3

\[ y_{n+1}=y_n+h y'_n+\dfrac{h^2}{2}y''_n+\dfrac{h^3}{6}y'''_n \]

where \(y'''=\dfrac{d}{dt}(f_t+f_y f)\) is built by repeated total differentiation. Order 3; stability function \(1+z+\dfrac{z^2}{2}+\dfrac{z^3}{6}\).

Taylor method, order 4 explicitorder 4

\[ y_{n+1}=y_n+h\,y'_n+\dfrac{h^2}{2}\,y''_n+\dfrac{h^3}{6}\,y'''_n+\dfrac{h^4}{24}\,y^{(4)}_n \]

The four total derivatives \(y',y'',y''',y^{(4)}\) are computed symbolically from \(f\). Order 4 with the same stability polynomial as RK4. (Requires \(f\) to be differentiable; otherwise this method is skipped with a note.)

Adams–Bashforth methods (explicit multistep)

Adams–Bashforth 2explicitorder 2

Explicit two-step method using the two most recent slopes.

\[ y_{n+1}=y_n+\dfrac h2\big(3f_n-f_{n-1}\big) \]

Order 2, one new \(f\) evaluation per step (the rest are reused). The single starting value \(y_1\) is generated with an explicit order-2 starter, improved Euler (Heun).

Adams–Bashforth 3 explicitorder 3

\[ y_{n+1}=y_n+\dfrac{h}{12}\big(23f_n-16f_{n-1}+5f_{n-2}\big) \]

Order 3. Explicit multistep regions of stability shrink with order, so AB methods need small \(h\) on stiff problems. The two starting values are generated with an explicit order-3 Runge–Kutta starter.

Adams–Bashforth 4 explicitorder 4

\[ y_{n+1}=y_n+\dfrac{h}{24}\big(55f_n-59f_{n-1}+37f_{n-2}-9f_{n-3}\big) \]

Order 4, four-step explicit method — efficient on smooth non-stiff problems since only one fresh \(f\) evaluation is needed per step. The three starting values are generated with classical RK4.

Adams–Moulton methods (implicit multistep)

Adams–Moulton 3implicitorder 3

Implicit two-step method (third order), solved with Newton each step.

\[ y_{n+1}=y_n+\dfrac{h}{12}\big(5f_{n+1}+8f_n-f_{n-1}\big) \]

For the same order, AM methods have much larger stability regions and smaller error constants than AB methods, at the cost of an implicit solve. The starting value is generated with an explicit order-3 Runge–Kutta starter.

Adams–Moulton 4 implicitorder 4

\[ y_{n+1}=y_n+\dfrac{h}{24}\big(9f_{n+1}+19f_n-5f_{n-1}+f_{n-2}\big) \]

Order 4 implicit three-step method. Often paired with AB4 as a predictor–corrector; here it is solved directly by Newton iteration. The two starting values are generated with classical RK4.

Adams–Moulton 5 implicitorder 5

\[ y_{n+1}=y_n+\dfrac{h}{720}\big(251f_{n+1}+646f_n-264f_{n-1}+106f_{n-2}-19f_{n-3}\big) \]

Order 5 implicit four-step method — the next Adams–Moulton formula. For a given order it has a much smaller error constant and a larger stability region than the explicit Adams–Bashforth method of the same order, at the cost of an implicit solve (done here by Newton iteration). Its three starting values are generated with an order-5 Runge–Kutta starter (Dormand–Prince) so the multistep method keeps its full order.

Backward differentiation formulas (BDF)

BDF 2 (backward differentiation)implicitorder 2

Backward differentiation formulas fit a polynomial through the new and past solution values and match its derivative to \(f\) at the new point — so only \(f(t_{n+1},y_{n+1})\) appears. They are the standard methods for stiff problems.

\[ y_{n+2}=\dfrac{4}{3}y_{n+1}-\dfrac{1}{3}y_n+\dfrac{2}{3}h\,f(t_{n+2},y_{n+2}) \]

Order 2 and A-stable — the only BDF that is. Solved by Newton each step. Its single starting value \(y_1\) is generated with Radau IIA (order 5, L-stable), an implicit starter that stays stable on stiff problems — matching BDF2's own intent.

BDF 3 implicitorder 3

\[ y_{n+3}=\dfrac{18}{11}y_{n+2}-\dfrac{9}{11}y_{n+1}+\dfrac{2}{11}y_n+\dfrac{6}{11}h\,f(t_{n+3},y_{n+3}) \]

Order 3. Not A-stable but A(\(\alpha\))-stable with \(\alpha\approx 86^\circ\): stable in a wide wedge of the left half-plane, which is enough for most stiff problems.

BDF 4 implicitorder 4

\[ y_{n+4}=\dfrac{48}{25}y_{n+3}-\dfrac{36}{25}y_{n+2}+\dfrac{16}{25}y_{n+1}-\dfrac{3}{25}y_n+\dfrac{12}{25}h\,f(t_{n+4},y_{n+4}) \]

Order 4, A(\(\alpha\))-stable with \(\alpha\approx 73^\circ\). The stability wedge narrows as the order rises.

BDF 5 implicitorder 5

\[ y_{n+5}=\dfrac{300}{137}y_{n+4}-\dfrac{300}{137}y_{n+3}+\dfrac{200}{137}y_{n+2}-\dfrac{75}{137}y_{n+1}+\dfrac{12}{137}y_n+\dfrac{60}{137}h\,f(t_{n+5},y_{n+5}) \]

Order 5, A(\(\alpha\))-stable with \(\alpha\approx 51^\circ\) — the highest-order BDF still usable (BDF 6 is unstable, BDF \(\ge 7\) are not zero-stable). Its four starting values are generated with Radau IIA (order 5, L-stable), whose order matches BDF5 exactly and stays stable on stiff problems.

Predictor–corrector methods

Adams PC 2 (ABM2)explicitorder 2

A predictor–corrector in PECE form: the explicit Adams–Bashforth 2 predictor supplies \(y^{*}_{n+1}\), it is evaluated once, and the implicit Adams–Moulton 2 (trapezoidal) corrector is applied a single time using that value — no implicit solve.

\[ \begin{aligned} y^{*}_{n+1}&=y_n+\dfrac{h}{2}\big(3f_n-f_{n-1}\big)\\ y_{n+1}&=y_n+\dfrac{h}{2}\big(f(t_{n+1},y^{*}_{n+1})+f_n\big) \end{aligned} \]

Order 2, two \(f\) evaluations per step. The single starting value is generated with improved Euler (Heun).

Adams PC 3 (ABM3) explicitorder 3

AB3 predictor with the AM3 corrector applied once (PECE).

\[ \begin{aligned} y^{*}_{n+1}&=y_n+\dfrac{h}{12}\big(23f_n-16f_{n-1}+5f_{n-2}\big)\\ y_{n+1}&=y_n+\dfrac{h}{12}\big(5f(t_{n+1},y^{*}_{n+1})+8f_n-f_{n-1}\big) \end{aligned} \]

Order 3. The two starting values come from a third-order Runge–Kutta starter.

Adams PC 4 (ABM4) explicitorder 4

The classic Adams–Bashforth–Moulton predictor–corrector: AB4 predicts, AM4 corrects once (PECE). It pairs the AB4 and AM4 methods above but, unlike standalone AM4, needs no Newton iteration — only two \(f\) evaluations per step.

\[ \begin{aligned} y^{*}_{n+1}&=y_n+\dfrac{h}{24}\big(55f_n-59f_{n-1}+37f_{n-2}-9f_{n-3}\big)\\ y_{n+1}&=y_n+\dfrac{h}{24}\big(9f(t_{n+1},y^{*}_{n+1})+19f_n-5f_{n-1}+f_{n-2}\big) \end{aligned} \]

Order 4 with a much larger stability region and smaller error constant than AB4 alone. Three starting values are generated with classical RK4.

Milne–Simpson explicitorder 4

Milne's fourth-order predictor (an open Newton–Cotes formula) with the Simpson-rule corrector.

\[ \begin{aligned} y^{*}_{n+1}&=y_{n-3}+\dfrac{4h}{3}\big(2f_n-f_{n-1}+2f_{n-2}\big)\\ y_{n+1}&=y_{n-1}+\dfrac{h}{3}\big(f(t_{n+1},y^{*}_{n+1})+4f_n+f_{n-1}\big) \end{aligned} \]

Order 4 with a tiny error constant, but weakly unstable: its characteristic polynomial has a parasitic root on the unit circle (\(w=-1\) at \(z=0\)), so errors fail to decay and the solution eventually oscillates — visible as an almost non-existent stability region. Three RK4 starting values.

Hamming explicitorder 4

Uses the same Milne predictor but Hamming's corrector, whose coefficients are chosen to push the parasitic roots strictly inside the unit circle — curing Milne–Simpson's weak instability.

\[ \begin{aligned} y^{*}_{n+1}&=y_{n-3}+\dfrac{4h}{3}\big(2f_n-f_{n-1}+2f_{n-2}\big)\\ y_{n+1}&=\dfrac{9y_n-y_{n-2}}{8}+\dfrac{3h}{8}\big(f(t_{n+1},y^{*}_{n+1})+2f_n-f_{n-1}\big) \end{aligned} \]

Order 4 and a stable, practical fourth-order multistep predictor–corrector. Three RK4 starting values.

Examples

To solve your own equation, choose Custom at the bottom of this list, then edit \(f(t,y)\), \(t_0\), \(y_0\) and \(t_{\text{end}}\) below.

Right-hand side \(f(t,y)\)

Start \(t_0\)

Initial \(y_0\)

End \(t_{\text{end}}\)

Parameter \(\lambda\)

−2.0

exact \(\lambda\):

Drag the slider or type an exact \(\lambda\) (any value); large negative makes \(y'=\lambda y\) stiff.

Adaptive step size

Number of steps \(N\)

N = 8

exact \(N\): h = 0.25

Live update (immediate)

Slide right to add steps (smaller \(h\)); each notch is one more step.

Graph height

380px

Global error norm:

Select methods and press Solve.

Solution & approximations

Export solution

About this plot

The equation and reference shown above the plot identify what is being solved. The black curve is the exact solution (for built-in examples) or a fine numerical reference (dense RK4, or an A-stable implicit solver when the problem is stiff). Colored curves with markers are the selected methods sampled at their step nodes. The view is auto-scaled to the reference so it fills the frame; divergent (unstable) runs simply shoot off the top or bottom — scroll to zoom or drag to pan to follow them.

About this plot

Absolute error at each grid point on a logarithmic axis. Higher-order methods sit lower; watch the error explode once a step leaves a method’s stability region. For the discontinuous right-hand sides (the staircase, relay, and square-wave examples) the error jumps at each discontinuity of the equation, so the line is broken there: it stays connected within each smooth stretch but leaves a gap across the step that straddles a jump, rather than drawing a misleading diagonal through it.

Computing…

Method:

About this plot

Shaded blue: \(z=h\lambda\) values for which the method is absolutely stable (\(|\)amplification\(|\le1\)); the dark curve is the stability boundary. The red × marks \(z=h\lambda\) for this problem. When no explicit \(\lambda\) is set, \(\lambda=\partial f/\partial y\) is sampled along the solution and the \(\times\) is placed at its most negative value — the worst case for absolute stability (hover to see at which \(t\)). If the × lies outside the shaded region, expect the method to grow there.

Computing…

About this plot

Log–log plot of the global error (in the norm chosen by the norm selector) against \(h\) (each \(h\) halved). The fitted empirical order \(p\) — the slope of the line — appears in the legend; a straight line of slope \(p\) confirms \(O(h^p)\) convergence. Step sizes where a method is unstable (large \(h\), error far beyond the solution scale) or round-off-dominated (very small \(h\), error stuck near machine precision) are dropped automatically, so the axis and the fitted order use only the clean convergent range. In adaptive mode this becomes a work–precision diagram: the error is plotted against the number of accepted steps \(N\) over a sweep of tolerances, and since error \(\approx C\,N^{-p}\) the slope of \(\log|\text{Error}|\) vs. \(\log N\) gives \(-p\); the fitted \(p\) shown in the legend is measured on the most-refined (tightest-tolerance) points.

Cite this tool

DOI: https://doi.org/10.5281/zenodo.20981245

Kapita, S. (2026). Initial Value Problem Methods: A Computational Laboratory. Math Tools. https://doi.org/10.5281/zenodo.20981245

Kapita, Shelvean. "Initial Value Problem Methods: A Computational Laboratory." Math Tools, 2026, doi.org/10.5281/zenodo.20981245.

@online{kapita2026ivp,
  author       = {Shelvean Kapita},
  title        = {{Initial Value Problem Methods: A Computational Laboratory}},
  year         = {2026},
  organization = {Math Tools},
  doi          = {10.5281/zenodo.20981245},
  url          = {https://doi.org/10.5281/zenodo.20981245}
}