Numerical Methods for ODEs

The local truncation error of a one-step method $y_{n+1} = y_n + h \cdot \Phi(t_n, y_n, h)$ at step $n$ is the residual: $\tau_n \;=\; y(t_{n+1}) - y(t_n) - h \cdot \Phi(t_n, y(t_n), h),$ where $y(t)$ is the exact solution of the IVP. This measures how closely a single step of the method matches the true Taylor expansion of $y$ — independent of any errors that have already accumulated in $y_n$.

Let $f: [t_0, T] \times \mathbb{R} \to \mathbb{R}$ be continuous and Lipschitz in $y$ with constant $L$, and assume the exact solution $y(t)$ is twice continuously differentiable with $|y''(t)| \leq M$ on $[t_0, T]$. Let $y_0, y_1, \ldots, y_N$ be the Euler approximations with step size $h$, and let $e_n = y_n - y(t_n)$ be the error at step $n$ (so $e_0 = 0$). Then $|e_n| \;\leq\; \frac{Mh}{2L}\bigl(e^{L(t_n - t_0)} - 1\bigr) \quad\text{for all } n.$ In particular, $|e_n| = O(h)$ as $h \to 0$, with the constant depending exponentially on $L$ and the integration interval.

We track how the error $e_n = y_n - y(t_n)$ evolves from one step to the next. By the definitions of the Euler step and the local truncation error, $y_{n+1} = y_n + h f(t_n, y_n)$ $y(t_{n+1}) = y(t_n) + h f(t_n, y(t_n)) + \tau_n$ Subtracting, $e_{n+1} = e_n + h\bigl[f(t_n, y_n) - f(t_n, y(t_n))\bigr] - \tau_n.$ Apply the Lipschitz condition $|f(t_n, y_n) - f(t_n, y(t_n))| \leq L|y_n - y(t_n)| = L|e_n|$ and the bound $|\tau_n| \leq Mh^2/2$ from the Taylor expansion above: $|e_{n+1}| \;\leq\; |e_n| + hL|e_n| + \tfrac{Mh^2}{2} \;=\; (1 + hL)|e_n| + \tfrac{Mh^2}{2}.$

This is a linear recurrence in $|e_n|$ with growth factor $(1 + hL)$ and forcing term $Mh^2/2$. Solving by induction (or by the discrete Grönwall lemma), starting from $e_0 = 0$: $|e_n| \;\leq\; \tfrac{Mh^2}{2} \sum_{k=0}^{n-1} (1 + hL)^k \;=\; \tfrac{Mh^2}{2} \cdot \frac{(1 + hL)^n - 1}{hL} \;=\; \tfrac{Mh}{2L}\bigl[(1 + hL)^n - 1\bigr].$

The geometric sum used the formula $\sum_{k=0}^{n-1} r^k = (r^n - 1)/(r - 1)$ with $r = 1 + hL$.

Now use the elementary inequality $1 + x \leq e^x$ for all real $x$, applied to $x = hL$: $(1 + hL)^n \;\leq\; (e^{hL})^n \;=\; e^{nhL} \;=\; e^{L(t_n - t_0)},$ since $nh = t_n - t_0$. Substituting, $|e_n| \;\leq\; \frac{Mh}{2L}\bigl(e^{L(t_n - t_0)} - 1\bigr).$

This bound has three pieces. The factor $Mh/(2L)$ shows the linear-in-$h$ scaling that defines first-order convergence. The factor $(e^{L(t_n - t_0)} - 1)$ is the error amplification — it grows exponentially in both the Lipschitz constant $L$ and the integration time. And the dependence on $M = \max|y''|$ shows that highly curved solutions accumulate error faster than nearly-linear ones, since the local truncation $\tau_n = (h^2/2)y''(\xi_n)$ depends on the second derivative. $\blacksquare$

Take $y' = y$, $y(0) = 1$, with exact solution $y(t) = e^t$. Here $f(t, y) = y$, so $L = 1$ and $y'' = y$, giving $M = e^T$ on $[0, T]$. The Euler iteration is $y_{n+1} = y_n + h y_n = (1+h) y_n$, so $y_n = (1+h)^n$.

At $T = 1$ with $h = 0.1$: $y_{10} = (1.1)^{10} \approx 2.594$, while the true value is $e \approx 2.718$. The error is $\approx 0.124$, which matches the bound $|e_n| \leq (Mh/2L)(e^L - 1) = (e \cdot 0.1 / 2)(e - 1) \approx 0.234$ (the bound is loose by a factor of 2 — Grönwall is conservative).

Halving the step to $h = 0.05$: $y_{20} = (1.05)^{20} \approx 2.653$, error $\approx 0.065$ — exactly half. This is the linear-in-$h$ scaling that defines a first-order method. To get an error of $10^{-6}$, you’d need $h \approx 10^{-6}$ — one million steps on the unit interval. Euler’s method is honest about its accuracy, but it makes you pay for it.

Now try $y' = -10y$, $y(0) = 1$, with exact solution $y(t) = e^{-10t}$. This is a fast-decaying equation: by $t = 1$, the true solution is $e^{-10} \approx 4.5 \times 10^{-5}$.

The Euler iteration is $y_{n+1} = (1 - 10h) y_n$, so $y_n = (1 - 10h)^n$. The behavior depends critically on $|1 - 10h|$:

• If $h = 0.01$: $1 - 10h = 0.9$, so $y_{100} = 0.9^{100} \approx 2.7 \times 10^{-5}$. Reasonably close to the true answer.
• If $h = 0.1$: $1 - 10h = 0$, so $y_n = 0$ for all $n \geq 1$. Wildly wrong (off by a factor of $\infty$ relatively, but at least bounded).
• If $h = 0.2$: $1 - 10h = -1$, so $y_n$ oscillates between $1$ and $-1$ forever. The numerical solution doesn’t even decay.
• If $h = 0.3$: $1 - 10h = -2$, so $y_n = (-2)^n$ — the numerical solution blows up, even though the true solution is decaying to zero!

This is shocking: making the step size larger didn’t just reduce accuracy — at $h = 0.3$, it caused the numerical solution to head to $\pm\infty$ while the true solution is going to zero. The problem isn’t truncation error: it’s stability. The Euler method amplifies errors by a factor of $|1 - 10h|$ each step, and this amplification factor exceeds 1 once $h > 0.2$. We’ll formalize this as the “stability region” of Euler’s method in Section 5.

Examples 1 and 2 reveal the fundamental tension in numerical methods. Truncation error pushes us to make $h$ small: Euler is first-order, so halving $h$ halves the error. But stability also constrains $h$: if $h$ is too large, errors in past steps get amplified into runaway oscillations or blow-up. For Euler on the test equation $y' = \lambda y$ with $\lambda < 0$, stability requires $h \leq 2/|\lambda|$.

For “well-behaved” equations these constraints don’t conflict — accuracy demands a smaller $h$ than stability does, and you just pick the smaller of the two. But for stiff equations, where $|\lambda|$ is large but $y(t)$ is also small, the stability constraint dominates and forces $h$ to be far smaller than accuracy requires. Section 6 develops this in detail. The cure is implicit methods, which have unconditional stability and break the tension. The cost: each step requires solving a nonlinear equation rather than evaluating one. Section 7 develops them.

A Runge-Kutta method is uniquely specified by its coefficients, conventionally arranged in a Butcher tableau (after John Butcher, who systematized the theory in the 1960s). The tableau has three pieces: a column vector $\mathbf c$ listing the time fractions of each stage, an $s \times s$ matrix $A$ holding the stage-coupling coefficients, and a row vector $\mathbf b$ holding the weights used to combine the stages into the final update. The conventional layout writes $\mathbf c$ to the left of $A$, with a horizontal rule separating the stage rows from the row of weights $\mathbf b$. For an explicit method, the matrix $A$ is strictly lower triangular ($a_{i,j} = 0$ for $j \geq i$), which means the stages can be computed in order $k_1, k_2, \ldots, k_s$ without any nonlinear solve.

Forward Euler corresponds to the trivial case $s = 1$, $c_1 = 0$, $a_{1,1} = 0$, $b_1 = 1$. The midpoint method uses $s = 2$ stages with $c_1 = 0$, $c_2 = 1/2$, $a_{2,1} = 1/2$, and weights $b_1 = 0$, $b_2 = 1$ to achieve order 2. The classical RK4 uses $s = 4$ stages with $c = (0, 1/2, 1/2, 1)$, with the matrix $A$ filled in so that $k_2$ uses $k_1$ at the half-step, $k_3$ uses $k_2$ at the half-step, and $k_4$ uses $k_3$ at the full step. The weights are $b_1 = 1/6$, $b_2 = 1/3$, $b_3 = 1/3$, $b_4 = 1/6$ — Simpson’s rule applied to the four stage-slopes.

The classical Runge-Kutta method defined above has local truncation error $O(h^5)$ and global error $O(h^4)$. Equivalently, it satisfies all eight of the order conditions required for a 4-stage explicit Runge-Kutta method to achieve order 4 in the Butcher theory.

Let’s compute one step of RK4 from $t_0 = 0$ with $h = 0.1$. The right-hand side is $f(t, y) = y$.

$k_1 = f(0, 1) = 1,$ $k_2 = f(0.05,\; 1 + 0.05 \cdot 1) = f(0.05, 1.05) = 1.05,$ $k_3 = f(0.05,\; 1 + 0.05 \cdot 1.05) = f(0.05, 1.0525) = 1.0525,$ $k_4 = f(0.1,\; 1 + 0.1 \cdot 1.0525) = f(0.1, 1.10525) = 1.10525,$ $y_1 = 1 + (0.1/6)(1 + 2(1.05) + 2(1.0525) + 1.10525) = 1 + 0.10517 = 1.10517.$

The exact value is $e^{0.1} = 1.10517091808\ldots$, so RK4 nailed it to 5 decimal places in a single step. Compare to Euler’s $y_1 = 1 + 0.1 \cdot 1 = 1.1$, which is correct only to 2 decimal places. RK4 used four function evaluations to get five extra digits — a staggering improvement per evaluation.

Suppose we want global error $\leq 10^{-3}$ on the IVP $y' = y$, $y(0) = 1$, integrated to $t = 1$. From the analysis above, Euler’s global error is $\sim h$ (with constant of order 1 here), so we need $h \approx 10^{-3}$, i.e., 1000 steps.

RK4’s global error is $\sim h^4$, so we need $h^4 \approx 10^{-3}$, i.e., $h \approx 0.18$. That’s about 6 steps. Each RK4 step costs 4 function evaluations, so RK4 uses $\approx 24$ function evaluations total. Euler uses $\approx 1000$. RK4 is 40 times faster for the same accuracy on this problem.

The crossover where Euler and RK4 cost the same per accuracy is roughly $h \approx 1$ — so for any reasonable accuracy target, RK4 wins. The story changes for stiff equations (Section 6) where stability, not accuracy, sets the step size; then the per-step cost of RK4 becomes a liability rather than an asset.

RK4 is the workhorse of explicit methods for moderate accuracy problems. Going higher — RK5, RK8 — gives diminishing returns. The number of order conditions grows combinatorially, so an order-8 method might need 11 stages to satisfy them all, costing 11 function evaluations per step. This pays off only if the savings in step count exceed the per-step overhead. For very high accuracy ($\varepsilon < 10^{-12}$) on smooth problems, methods like Dormand-Prince 8(7) or extrapolation methods become competitive. For low accuracy ($\varepsilon > 10^{-3}$), Euler or RK2 may suffice. RK4 hits the sweet spot for typical problems.

The right way to think about it: the “order vs. cost” trade-off depends on (a) the accuracy you need, (b) the cost of one $f$ evaluation, and (c) the smoothness of $f$. For neural ODEs, $f$ is a forward pass through a network — expensive — and the desired accuracy is moderate, so RK4 or RK45 are both popular choices. For molecular dynamics with millions of atoms, $f$ is even more expensive, and methods specialized for Hamiltonian structure (symplectic integrators) take precedence over generic high-order methods.

An embedded RK pair of orders $p$ and $p+1$ consists of two sets of weights $\{b_i\}$ and $\{\hat b_i\}$ applied to the same $s$ stages $k_1, \ldots, k_s$: $y_{n+1} = y_n + h \sum_i b_i k_i \quad (\text{order } p), \qquad \hat y_{n+1} = y_n + h \sum_i \hat b_i k_i \quad (\text{order } p+1).$ The difference $\hat y_{n+1} - y_{n+1} = h \sum_i (\hat b_i - b_i) k_i$ is an estimate of the local truncation error of the order-$p$ method, accurate to leading order. The most widely used embedded pair is the Dormand-Prince RK4(5) method (Dormand & Prince, 1980), which has $s = 7$ stages and is the default solver in scipy.integrate.solve_ivp (method='RK45') and torchdiffeq ('dopri5').

For an embedded RK pair where the lower-order method has order $p$, the local truncation error scales like $\text{errEst} \approx C h^{p+1}$ for some constant $C$. To achieve a target tolerance $\text{tol}$ on the next step, we want $\text{tol} = C \tilde h^{p+1}$, where $\tilde h$ is the new step size. Dividing the two equations and solving: $\tilde h \;=\; h \cdot \left(\frac{\text{tol}}{\text{errEst}}\right)^{1/(p+1)}.$ In practice this is multiplied by a safety factor of about $0.9$ (so the new step is slightly conservative), and clamped to a range $[h_{\min}, h_{\max}]$ to prevent runaway shrinkage or growth. If $\text{errEst} > \text{tol}$, the current step is rejected (the computed $y_{n+1}$ is discarded) and the step is retried with the smaller $\tilde h$. If $\text{errEst} \leq \text{tol}$, the step is accepted and the larger $\tilde h$ is used for the next step.

Consider the forced linear ODE $y' = -50(y - \cos t) - \sin t, \qquad y(0) = 0,$ which has the closed-form solution $y(t) = \cos t \cdot (1 - e^{-50 t})$ — a sinusoid that “rings up” to its steady-state amplitude over a transient of length $1/50$. For $t \lesssim 0.1$ the solution changes rapidly (rising from 0 to nearly $\cos t$); for $t \gtrsim 0.1$ it follows the smooth cosine.

A fixed-step RK4 with $h = 0.01$ would take $100$ steps to reach $t = 1$, mostly in the smooth region where $h = 0.05$ would have sufficed. RK45 with $\text{tol} = 10^{-6}$ uses small steps in the transient ($h \sim 0.005$ for the first few accepted steps) and then expands rapidly: by $t = 0.5$ it’s using $h \sim 0.1$. Total step count: maybe $\sim 30$ — three times fewer than the fixed-step version, with comparable accuracy. The savings come from concentrating computational effort where the solution actually demands it.

The IVP $y' = y^2$, $y(0) = 1$ has the explicit solution $y(t) = 1/(1-t)$, which blows up as $t \to 1^-$. A fixed-step method either misses the blow-up entirely (if $h$ is too small to reach $t = 1$) or overshoots and produces nonsense.

RK45 handles this gracefully: as $y$ grows, $|f(t, y)| = y^2$ grows even faster, and the local truncation error explodes proportionally. The step-size formula $\tilde h = h \cdot (\text{tol}/\text{errEst})^{1/5}$ then shrinks $h$ rapidly. Near the singularity, RK45 takes geometrically smaller steps until $h$ hits the floor $h_{\min}$, at which point the integration halts. The user gets a reasonable solution up to $t \approx 0.99$ or so, with a clear signal (extreme step-size shrinkage, or hitting $h_{\min}$) that something is going wrong. This is the “graceful degradation” that adaptive methods provide and fixed-step methods don’t.

In practice, scientific ODE solvers use a combined tolerance: $\text{tol}_n = \text{atol} + \text{rtol} \cdot |y_n|$. The absolute tolerance $\text{atol}$ kicks in when $|y_n|$ is small (close to the noise floor), and the relative tolerance $\text{rtol}$ kicks in when $|y_n|$ is large (so the error is bounded as a fraction of the value). Typical defaults are $\text{atol} = 10^{-6}$ and $\text{rtol} = 10^{-3}$.

This dual scheme avoids two failure modes: a pure relative tolerance breaks down when the solution crosses zero ($\text{tol} \to 0$ would force $h \to 0$); a pure absolute tolerance is too coarse for small solutions and too tight for large ones. In scipy.integrate.solve_ivp, you pass atol and rtol as arguments. In torchdiffeq, the same conventions apply. ML practitioners using neural ODEs should be aware of these — the wrong choice of tolerance can either slow training to a crawl (over-tight) or introduce solver-induced bias into gradients (over-loose).

The stability region of a numerical method is the set $S \;=\; \{ z \in \mathbb{C} : |R(z)| \leq 1 \},$ where $R(z)$ is the method’s stability function obtained by applying the method to $y' = \lambda y$ with $z = h\lambda$. The numerical solution decays for an eigenvalue $\lambda$ if and only if $h\lambda \in S$.

For a scalar real $\lambda < 0$, this gives the constraint $h \leq h_{\max}(\lambda)$ where $h_{\max}$ is determined by where the boundary of $S$ crosses the real axis. For a system $y' = Ay$, every eigenvalue of $hA$ must lie in $S$, which is the constraint that determines the maximum stable step size for a multi-mode system.

A numerical method is A-stable if its stability region $S$ contains the entire closed left half-plane $\mathbb{C}^- = \{z : \text{Re}(z) \leq 0\}$. Equivalently, $|R(z)| \leq 1$ for all $z$ with $\text{Re}(z) \leq 0$.

A-stability is the strongest possible stability property: an A-stable method gives a decaying numerical solution for any stable scalar test equation, regardless of step size. There are no step-size restrictions for stability — only for accuracy.

There is no explicit Runge-Kutta method that is A-stable. More generally, no explicit linear multistep method of order greater than $2$ can be A-stable, and no implicit linear multistep method of order greater than $2$ can be A-stable. (The trapezoidal rule, which is implicit and second-order, is A-stable — it sits exactly at the barrier.)

Apply forward Euler to $y' = \lambda y$: $y_{n+1} = y_n + h \lambda y_n = (1 + h\lambda) y_n.$ So $R(z) = 1 + z$, where $z = h\lambda$. The stability region is $S = \{z \in \mathbb{C} : |1 + z| \leq 1\}.$ This is the closed disk of radius 1 centered at $z = -1$. On the real axis, it’s the interval $[-2, 0]$. So for a real eigenvalue $\lambda < 0$, forward Euler is stable iff $h \leq 2/|\lambda|$. This recovers the constraint from Example 2 in Section 2: on $y' = -10y$, stability requires $h \leq 0.2$ — exactly where we saw the numerical solution start oscillating.

Applying RK4 to $y' = \lambda y$ and computing the four stages, the stability function turns out to be exactly the truncated Taylor series of $e^z$ to fourth order: $R(z) = 1 + z + \tfrac{z^2}{2} + \tfrac{z^3}{6} + \tfrac{z^4}{24}.$ This is no accident: applying RK4 to the test equation forces $R$ to match $e^z = e^{h\lambda}$, the exact solution, to as many Taylor terms as the method’s order — for RK4 that’s 5 terms (orders 0 through 4).

The stability region $\{z : |R(z)| \leq 1\}$ is now a four-lobed shape that extends along the real axis to about $z \approx -2.78$. So RK4 is stable for real $\lambda < 0$ when $h < 2.78/|\lambda|$ — a roughly 40% improvement over forward Euler, in addition to RK4’s much higher accuracy. RK4 is not A-stable: its region is bounded, so for very stiff systems it still needs tiny step sizes for stability.

We compute the stability region directly. Forward Euler applied to $y' = \lambda y$ gives $y_{n+1} = (1 + h\lambda) y_n$, so iteratively $y_n = (1 + h\lambda)^n y_0 = R(z)^n y_0$ where $z = h\lambda$.

The numerical solution remains bounded iff $|R(z)| \leq 1$, which is iff $|1 + z| \leq 1$. Writing $z = x + iy$ for real $x, y$: $|1 + z|^2 = (1 + x)^2 + y^2 \leq 1 \iff (1 + x)^2 + y^2 \leq 1.$ This is the closed disk of radius $1$ centered at $z = -1$. It is a compact subset of $\mathbb{C}$, contained entirely in $\{z : -2 \leq \text{Re}(z) \leq 0\}$.

The disk does not contain the entire left half-plane: for example, $z = -10$ gives $|1 + z| = 9 > 1$. Therefore forward Euler is not A-stable. For any real eigenvalue $\lambda < 0$, the constraint $h\lambda \in S$ becomes $-2 \leq h\lambda \leq 0$, i.e., $h \leq 2/|\lambda|$. As $|\lambda|$ grows, the maximum stable step size shrinks proportionally — this is the source of stiffness pain for explicit methods. $\blacksquare$

For a system $\dot{\mathbf{y}} = \mathbf{f}(\mathbf{y})$, the stiffness ratio at a point $\mathbf{y}$ is $\kappa(\mathbf{y}) \;=\; \frac{\max_i |\lambda_i(J(\mathbf{y}))|}{\min_i |\lambda_i(J(\mathbf{y}))|},$ where $\lambda_i$ are the eigenvalues of the Jacobian $J(\mathbf{y}) = \partial \mathbf{f}/\partial \mathbf{y}$. A system is stiff at $\mathbf{y}$ if $\kappa(\mathbf{y}) \gg 1$ — typically $\kappa > 100$ or so. Stiffness depends on the integration interval: a system with $\kappa = 10^6$ is stiff for $T = 1$ (where you need to integrate through $\sim 10^6/T$ fast-decay periods) but not necessarily stiff for $T = 10^{-7}$ (where the slow dynamics haven’t even started).

Note the connection to Topic 23: the Jacobian eigenvalues are exactly what we computed for linearization stability analysis. A stable equilibrium has $\text{Re}(\lambda_i) < 0$ for all eigenvalues; a stiff stable equilibrium has some eigenvalues with very large $|\text{Re}(\lambda_i)|$ alongside others with small $|\text{Re}(\lambda_i)|$.

A canonical stiff test problem: three coupled chemical reactions $\dot y_1 = -0.04\, y_1 + 10^4\, y_2 y_3,$ $\dot y_2 = 0.04\, y_1 - 10^4\, y_2 y_3 - 3 \times 10^7\, y_2^2,$ $\dot y_3 = 3 \times 10^7\, y_2^2,$ with $\mathbf{y}(0) = (1, 0, 0)$. The three rate constants span 9 orders of magnitude: $0.04$, $10^4$, $3 \times 10^7$. The Jacobian eigenvalues at the initial state have moduli ranging from $0.04$ to $\sim 10^7$, giving a stiffness ratio $\kappa \sim 10^9$.

For an explicit method on this problem, stability requires $h \cdot 10^7 \lesssim 1$, i.e., $h \lesssim 10^{-7}$. To integrate to $t = 100$ (a typical time scale for the slow $y_1$ component), you’d need $\sim 10^9$ steps. With each step requiring 4 function evaluations for RK4, that’s $\sim 4 \times 10^9$ function evaluations — completely intractable. An implicit method (which we develop in Section 7) can use $h \sim 0.1$ throughout, requiring only $\sim 1000$ steps. The speedup is six orders of magnitude.

The Van der Pol oscillator $\ddot x - \mu(1 - x^2)\dot x + x = 0$, written as a first-order system with $y_1 = x$, $y_2 = \dot x$: $\dot y_1 = y_2, \qquad \dot y_2 = \mu(1 - y_1^2) y_2 - y_1.$ For small $\mu$ this is a smooth oscillator. For large $\mu$ — say $\mu = 1000$ — it becomes a relaxation oscillator: long stretches where the solution is nearly constant, punctuated by very rapid jumps when $|y_1|$ crosses $1$ and the damping coefficient $\mu(1 - y_1^2)$ flips sign.

During the slow stretches the Jacobian’s eigenvalues are small. During the rapid jumps, one eigenvalue becomes very large (roughly $-2\mu$ at the most stiff points). The stiffness ratio varies from $O(1)$ during smooth stretches to $O(\mu^2) = O(10^6)$ during jumps. An explicit method has to use a step size suited to the worst stretches throughout, even though most of the trajectory is smooth. Adaptive step-size control helps somewhat (it can take big steps during smooth stretches and shrink during jumps), but even RK45 spends most of its computation in the rapid-jump regions and is still much slower than an implicit method. The Van der Pol oscillator at $\mu = 1000$ is a textbook case where you should reach for solve_ivp(method='Radau') or 'BDF' rather than the default 'RK45'.

A subtle but important point: stiffness is a combination of the equation, the initial condition, and the desired integration interval. A system with $\kappa = 10^6$ is stiff only if you want to integrate through many fast-decay periods. If your interval is short enough that the fast modes are still active, you actually need to resolve them — and an explicit method with small $h$ is appropriate.

This is why “stiff” is not a binary classification but a regime. Most modern ODE solvers (solve_ivp, DifferentialEquations.jl) include automatic stiffness detection: they monitor the Jacobian eigenvalues (or proxies like the step-size rejection rate of an explicit method) and switch between explicit and implicit methods on the fly. The field of ODE solving has matured to the point where users mostly don’t need to make this choice manually — but understanding why stiffness matters lets you debug when the automation gets it wrong.

Implicit (backward) Euler: $y_{n+1} = y_n + h \cdot f(t_{n+1}, y_{n+1})$. The unknown $y_{n+1}$ appears inside $f$ on the right-hand side — this is an implicit equation that must be solved at each step. For linear $f$, the solve is a matrix inversion; for nonlinear $f$, it’s a Newton iteration. This method is first-order accurate (same order as forward Euler) but A-stable.

Trapezoidal rule: $y_{n+1} = y_n + (h/2)[f(t_n, y_n) + f(t_{n+1}, y_{n+1})]$. The average of forward and backward Euler. Second-order accurate and A-stable — sitting exactly at the Dahlquist barrier of “highest-order A-stable linear multistep method.” Known as the Crank-Nicolson method in the PDE literature.

The implicit Euler method has stability function $R(z) = \frac{1}{1 - z},$ and is A-stable: $|R(z)| \leq 1$ for all $z \in \mathbb{C}$ with $\text{Re}(z) \leq 0$.

Apply implicit Euler to the test equation $y' = \lambda y$: $y_{n+1} = y_n + h\lambda y_{n+1}.$ Solving algebraically for $y_{n+1}$ (this is one of the rare cases where the implicit equation has a closed-form solution): $(1 - h\lambda) y_{n+1} = y_n \implies y_{n+1} = \frac{1}{1 - h\lambda} y_n = \frac{1}{1 - z} y_n,$ where $z = h\lambda$. So the stability function is $R(z) = 1/(1 - z)$.

To check A-stability, we need $|R(z)| \leq 1$ for all $z$ with $\text{Re}(z) \leq 0$. Write $z = x + iy$ with $x \leq 0$. Then $1 - z = (1 - x) - iy$, so $|1 - z|^2 = (1 - x)^2 + y^2.$ Since $x \leq 0$, we have $1 - x \geq 1$, so $(1 - x)^2 \geq 1$. Therefore $|1 - z|^2 = (1 - x)^2 + y^2 \geq 1 + y^2 \geq 1,$ which gives $|1 - z| \geq 1$ and so $|R(z)| = \frac{1}{|1 - z|} \leq 1.$ Equality $|R(z)| = 1$ holds iff $x = 0$ and $y = 0$ — that is, only at $z = 0$. Everywhere else in the closed left half-plane, $|R(z)| < 1$ strictly. So implicit Euler is A-stable, and in fact strongly damping for any $\text{Re}(\lambda) < 0$. $\blacksquare$

On the Robertson kinetics from Example 9, implicit Euler (with Newton iteration on the $3 \times 3$ Jacobian) can use $h = 10^{-2}$ to $h = 1$ throughout the integration. To reach $t = 100$, that’s $\sim 100$ to $10^4$ steps — six to seven orders of magnitude fewer than the $10^9$ steps an explicit method would need. Each implicit step is more expensive (Newton iteration costs $\sim 4 \times$ what one Euler step does), but $4 \cdot 10^4$ is still vastly less than $10^9$.