Chapter 23: Ordinary Differential Equations — Euler, Runge-Kutta, and Simulating Change

"The profound study of nature is the most fertile source of mathematical discoveries." — Joseph Fourier, The Analytical Theory of Heat (1822)

Overview

Almost everything that changes in time obeys a differential equation. A planet's orbit, the decay of a radioisotope, the concentration of a reactant, the population of a predator and its prey, the current in a circuit, the temperature of a cooling plate — each is governed by a rule that says how fast the state changes as a function of the state itself. Write that rule down and you have an ordinary differential equation. Solve it and you can predict the future: where the planet will be next year, how much of the isotope survives to next century, whether the population settles or oscillates. Most such equations have no closed-form solution, which is exactly why they belong in a book about scientific computing. We do not solve them with a pencil; we march them forward, one small step at a time, on a machine.

This chapter is about doing that marching correctly and efficiently in Fortran. We start with the simplest possible method — Euler's — precisely because it is simple enough to see through, and then we watch it fail: it is accurate only to first order, and buying accuracy by shrinking the step is ruinously expensive. The fix is the fourth-order Runge-Kutta method, RK4, the workhorse that has integrated more trajectories than any other algorithm in scientific history and the one you will reach for by default. From there we make the integration adaptive — letting the code choose its own step size — and vector-valued, so it handles systems of equations, which is what real problems always are.

Then comes the chapter's pivot, and the reason it sits where it does in the book. Discretize the space in a partial differential equation but leave time continuous, and the PDE collapses into an enormous system of ordinary differential equations — one equation per grid point. This is the method of lines, and it is the bridge from this chapter to the next. It also means that the heat-solver time-stepping you have been assembling since Chapter 4 is, viewed correctly, nothing but Euler's method applied to a very large ODE system — and that a better time integrator, like RK4, could be dropped straight in.

In this chapter, you will learn to:

• Set up an initial-value problem and write its right-hand side as a Fortran procedure you pass to a solver.
• Implement Euler's method, RK4, and an adaptive stepper, and reason about their accuracy by hand.
• Integrate systems of ODEs — orbits, oscillators, reaction networks — as vector states.
• See a PDE as a giant system of ODEs (the method of lines), and connect it to your heat solver.
• Recognize a stiff problem and understand why it demands an implicit method.

Learning Paths

How to read this chapter by track. - 🔬 Scientist — this is a chapter for you; read §23.1–23.4 closely and §23.6 for your domain. The method of lines in §23.4 is the idea that unifies your ODE and PDE work. - 📖 Standard — read straight through. Each method builds on the last, and the RK4 you meet in §23.2 is the single most useful algorithm in the chapter. - 🔧 Legacy — the classic FORTRAN 77 numerical libraries (ODEPACK, the rk routines of Numerical Recipes) implement exactly these methods; §23.2 and §23.5 explain what that inherited code is doing. - ⚡ HPC — skim §23.1–23.3; spend your time on §23.4 (systems and the method of lines) and §23.5 (stiffness), because they determine whether your time-stepping parallelizes well and how small your step must be.

23.1 The Initial-Value Problem and Euler's Method

An ordinary differential equation (ODE) relates a function to its own derivatives with respect to a single variable — here, time $t$. The problems we solve numerically almost always come packaged as an initial-value problem: an ODE together with the state at a starting time.

Definition (initial-value problem). A first-order initial-value problem (IVP) is the pair $$\frac{dy}{dt} = f(t, y), \qquad y(t_0) = y_0.$$ The function $f(t, y)$ — the right-hand-side function, or RHS — gives the instantaneous rate of change of the state $y$ as a function of the current time $t$ and current state $y$. The initial condition $y(t_0) = y_0$ pins down which particular solution we want out of the infinite family the ODE admits. To solve the IVP numerically is to produce approximate values $y_1, y_2, \dots$ of the true solution $y(t)$ at times $t_1, t_2, \dots$ marching forward from $t_0$.

The RHS is the whole physics of the problem. For radioactive decay, $f(t, y) = -\lambda y$ (the state decreases in proportion to how much is left). For our running test problem, $f(t, y) = y$ — the state grows in proportion to itself — whose exact solution is $y(t) = e^t$, a function we know to as many digits as we like and can therefore check our numerics against mercilessly.

The geometric picture is worth holding onto. At every point $(t, y)$ in the plane, $f$ hands you a slope — the direction the solution is heading. The exact solution is the curve that is everywhere tangent to that field of slopes and passes through $(t_0, y_0)$. Every method in this chapter is a different strategy for following that field forward without knowing the curve in advance.

Euler's method: follow the slope you have

The most obvious strategy is the right place to start. You are standing at $(t_n, y_n)$. The ODE tells you the slope there is $f(t_n, y_n)$. So take a small step $h$ forward in time and follow that slope in a straight line:

Definition (Euler's method). Euler's method advances the solution by $$y_{n+1} = y_n + h\, f(t_n, y_n), \qquad t_{n+1} = t_n + h,$$ where $h$ is the step size. It approximates the solution over each interval by the tangent line at the interval's left endpoint. It is the simplest ODE integrator that exists, and every more accurate method can be understood as a correction to it.

Where does it come from? Taylor's theorem: $y(t_n + h) = y(t_n) + h\,y'(t_n) + \tfrac{1}{2}h^2 y''(\xi)$. Euler keeps the first two terms — and $y'(t_n)$ is exactly $f(t_n, y_n)$ — and throws away the rest. The piece it throws away, $\tfrac{1}{2}h^2 y''(\xi)$, is the local truncation error of a single step: it is $O(h^2)$. (This is the truncation error you met for finite differences in Chapter 22, now applied to stepping forward in time instead of estimating a derivative.)

Here is Euler's method in Fortran. Notice how we pass the RHS as a procedure argument, using an abstract interface to describe its shape — the same procedure machinery from Chapter 6, now letting one solver work for any equation you hand it:

! example-01-euler.f90 -- Euler's method for the IVP y' = y, y(0) = 1
module euler_mod
  implicit none
  integer, parameter :: dp = selected_real_kind(15, 307)

  abstract interface
    function rhs(t, y) result(dydt)   ! the shape every RHS must have
      import :: dp
      real(dp), intent(in) :: t, y
      real(dp)             :: dydt
    end function rhs
  end interface

contains

  function euler_step(f, t, y, h) result(y_next)
    procedure(rhs)       :: f          ! a dummy procedure: the RHS f(t, y)
    real(dp), intent(in) :: t, y, h
    real(dp)             :: y_next
    y_next = y + h * f(t, y)
  end function euler_step

  function grow(t, y) result(dydt)     ! our test RHS: f(t, y) = y
    real(dp), intent(in) :: t, y
    real(dp)             :: dydt
    dydt = y
  end function grow
end module euler_mod

program run_euler
  use euler_mod
  implicit none
  real(dp) :: t, y, h
  integer  :: n

  h = 0.25_dp;  t = 0.0_dp;  y = 1.0_dp
  print '(a)', '   n           t            y'
  print '(i4, 2f13.6)', 0, t, y
  do n = 1, 4
    y = euler_step(grow, t, y, h)
    t = t + h
    print '(i4, 2f13.6)', n, t, y
  end do
  print '(a, f13.6)', 'exact  e = ', exp(1.0_dp)
  print '(a, f13.6)', 'error    = ', exp(1.0_dp) - y
end program run_euler

$ gfortran -std=f2018 -Wall example-01-euler.f90 -o euler && ./euler
   n           t            y
   0     0.000000     1.000000
   1     0.250000     1.250000
   2     0.500000     1.562500
   3     0.750000     1.953125
   4     1.000000     2.441406
exact  e =      2.718282
error    =      0.276876

Trace it by hand and you will see exactly what Euler does. Each step multiplies $y$ by $(1 + h) = 1.25$, because $y_{n+1} = y_n + h y_n = (1+h)y_n$. So after four steps $y = 1.25^4 = 2.44140625$ — and the true answer is $e = 2.71828\ldots$. Euler undershoots by more than ten percent, because at every step it follows the slope at the left end of the interval, and for a curve that is bending upward that slope is always too shallow. The straight-line tangent falls below the curve, and the error compounds.

💡 Intuition: Euler is a driver who checks the road direction, then closes their eyes for the whole step and steers straight. If the road is straight (the solution is linear) they are fine. If it curves, they drift to the outside of every bend, and the drift accumulates. Every better method in this chapter is a way of peeking at the road during the step.

The cost of first-order accuracy

The local error of one step is $O(h^2)$, but that is not the number that matters. To integrate from $t_0$ to a fixed final time $T$ you take $N = (T - t_0)/h$ steps, and the local errors accumulate. Very roughly, $N$ steps of $O(h^2)$ error each give a total, or global, error of $N \cdot O(h^2) = \tfrac{T - t_0}{h}\cdot O(h^2) = O(h)$.

Definition (order of accuracy for an integrator). A time integrator has global order $p$ if its error at a fixed final time shrinks like $O(h^p)$ as $h \to 0$. Euler's method is first order ($p = 1$): halve the step and you roughly halve the error. This is the same notion of order of accuracy you measured for quadrature in Chapter 22 — now it governs how fast a simulation converges as you refine the timestep.

First order is a hard bargain. Watch what it costs. Integrating $y' = y$ to $t = 1$ and comparing the error at several step sizes:

Every halving of $h$ roughly halves the error — the error ratios march toward $2$ — which is the signature of first order. To get one more correct decimal digit (a $10\times$ error reduction) you need about $10\times$ as many steps, and each step is a full evaluation of $f$. For a smooth problem this is an appalling deal, and it is the entire motivation for the next section.

🐍 Python Comparison: In Python you would never write Euler by hand for real work; you would call scipy.integrate.solve_ivp(f, [0, 1], [1.0]), which defaults to an adaptive fifth-order Runge-Kutta. But solve_ivp calls your Python f back on every stage, and if f is itself Python, that callback is the bottleneck — this is precisely the case where a Fortran RHS, wrapped for Python as in Chapter 15, pays off. The integrator's structure is identical; the speed is not.

🔄 Check Your Understanding. 1. In the IVP $y' = f(t, y),\ y(t_0) = y_0$, what does $f$ physically represent, and what does the initial condition do? 2. Euler's local truncation error is $O(h^2)$ but its global error is only $O(h)$. Where did the extra power of $h$ go? 3. You halve the step size in an Euler integration to a fixed final time. Roughly what happens to the error, and roughly what happens to the run time?

Answers

1. $f$ is the instantaneous rate of change (the slope field); the initial condition selects the one solution curve passing through $(t_0, y_0)$. 2. Reaching a fixed $T$ takes $N = (T-t_0)/h \propto 1/h$ steps, so the $N$ accumulated local errors of size $O(h^2)$ sum to $O(h)$. 3. The error roughly halves (first order); the run time roughly doubles (twice as many steps). One more correct digit costs about ten times the work.

23.2 Runge-Kutta: RK4, the Workhorse

Euler's flaw is that it commits to a single slope — the one at the start of the step — and lives with it. The Runge-Kutta idea, from Carl Runge and Wilhelm Kutta around 1900, is to sample the slope at several points within the step and combine those samples into a much better average slope. The samples cost extra evaluations of $f$, but each one buys a higher power of $h$ in accuracy, and that trade is overwhelmingly worth it.

Definition (Runge-Kutta method). A Runge-Kutta method advances one step by evaluating the RHS at a set of intermediate stages within $[t_n, t_n + h]$ and taking a weighted average of those stage slopes as the effective slope for the step. Different choices of stages and weights give methods of different order. The idea is to match more terms of the Taylor expansion of the true solution than Euler's single term.

Start with the simplest improvement, the midpoint method (a second-order, or RK2, method): use Euler to peek at the middle of the step, read the slope there, and use that slope for the whole step.

$$ k_1 = f(t_n, y_n), \qquad k_2 = f\!\left(t_n + \tfrac{h}{2},\ y_n + \tfrac{h}{2}k_1\right), \qquad y_{n+1} = y_n + h\,k_2. $$

On $y' = y$ with $h = 1$: $k_1 = 1$, then $k_2 = f(0.5,\ 1 + 0.5) = 1.5$, so $y_1 = 1 + 1(1.5) = 2.5$. Euler gave $2.0$; the midpoint method gives $2.5$; the true value is $2.718\ldots$. One extra evaluation of $f$ closed most of the gap. That is the Runge-Kutta bargain in miniature.

The classical fourth-order method

The method you will actually use samples four stages and is fourth-order accurate. It is the RK4, the classical Runge-Kutta method, and its coefficients are worth committing to memory because you will type them for the rest of your career:

$$ \begin{aligned} k_1 &= f\!\left(t_n,\ y_n\right) \\ k_2 &= f\!\left(t_n + \tfrac{h}{2},\ y_n + \tfrac{h}{2}\,k_1\right) \\ k_3 &= f\!\left(t_n + \tfrac{h}{2},\ y_n + \tfrac{h}{2}\,k_2\right) \\ k_4 &= f\!\left(t_n + h,\ y_n + h\,k_3\right) \\ y_{n+1} &= y_n + \frac{h}{6}\left(k_1 + 2k_2 + 2k_3 + k_4\right) \end{aligned} $$

Read the structure. Stage $k_1$ is the slope at the left end (Euler's slope). Stages $k_2$ and $k_3$ are two successive estimates of the slope at the midpoint, each using the previous stage to get there. Stage $k_4$ is the slope at the right end, reached using $k_3$. The final combination is a weighted average that trusts the two midpoint slopes twice as much as the endpoints — weights $\tfrac{1}{6}, \tfrac{2}{6}, \tfrac{2}{6}, \tfrac{1}{6}$, which sum to one, exactly the pattern of Simpson's rule for integration. That is not a coincidence: integrating $y' = f$ over the step is quadrature, and RK4's weights are Simpson's weights.

The compact bookkeeping for any Runge-Kutta method is the Butcher tableau — the stage times down the left, the coupling coefficients in the body, the final weights along the bottom. For RK4 it is:

   0  |
  1/2 | 1/2
  1/2 |  0    1/2
   1  |  0     0     1
------+-------------------------
      | 1/6   1/3   1/3   1/6

⚠️ Common Pitfall — the coefficients are exact, and they are unforgiving. The single most common RK4 bug is a mistyped coefficient: a $k_2$ that reuses $k_1$ where it should use $k_2$, a final weight of $\tfrac{1}{3}$ where it should be $\tfrac{1}{6}$, or an $h/2$ dropped to $h$. The program still compiles, still runs, still produces plausible-looking numbers — and is silently second-order or worse. Type the four stages exactly as above, and test the order numerically (halve $h$, confirm the error drops by about $16$); a wrong coefficient shows up immediately as the wrong convergence rate.

Here is RK4 in Fortran, reusing the procedure-argument pattern so the stepper is independent of the equation. We take one step of size $h = 1$ so every stage is a round number you can check against the hand trace:

! example-02-rk4.f90 -- classical RK4 for y' = y, one step of h = 1 from y(0) = 1
module rk4_mod
  implicit none
  integer, parameter :: dp = selected_real_kind(15, 307)

  abstract interface
    function rhs(t, y) result(dydt)
      import :: dp
      real(dp), intent(in) :: t, y
      real(dp)             :: dydt
    end function rhs
  end interface

contains

  function rk4_step(f, t, y, h, k1, k2, k3, k4) result(y_next)
    procedure(rhs)        :: f
    real(dp), intent(in)  :: t, y, h
    real(dp), intent(out) :: k1, k2, k3, k4   ! returned so we can print them
    real(dp)              :: y_next
    k1 = f(t,             y)
    k2 = f(t + 0.5_dp*h,  y + 0.5_dp*h*k1)
    k3 = f(t + 0.5_dp*h,  y + 0.5_dp*h*k2)
    k4 = f(t + h,         y + h*k3)
    y_next = y + (h/6.0_dp) * (k1 + 2.0_dp*k2 + 2.0_dp*k3 + k4)
  end function rk4_step

  function grow(t, y) result(dydt)
    real(dp), intent(in) :: t, y
    real(dp)             :: dydt
    dydt = y
  end function grow
end module rk4_mod

program run_rk4
  use rk4_mod
  implicit none
  real(dp) :: y1, k1, k2, k3, k4

  y1 = rk4_step(grow, 0.0_dp, 1.0_dp, 1.0_dp, k1, k2, k3, k4)
  print '(a, f12.6)', 'k1        = ', k1
  print '(a, f12.6)', 'k2        = ', k2
  print '(a, f12.6)', 'k3        = ', k3
  print '(a, f12.6)', 'k4        = ', k4
  print '(a, f12.6)', 'RK4  y(1) = ', y1
  print '(a, f12.6)', 'exact  e  = ', exp(1.0_dp)
  print '(a, f12.6)', 'RK4 error = ', exp(1.0_dp) - y1
end program run_rk4

$ gfortran -std=f2018 -Wall example-02-rk4.f90 -o rk4 && ./rk4
k1        =     1.000000
k2        =     1.500000
k3        =     1.750000
k4        =     2.750000
RK4  y(1) =     2.708333
exact  e  =     2.718282
RK4 error =     0.009948

Follow the stages: $k_1 = 1$; $k_2 = f(0.5,\ 1 + 0.5\cdot1) = 1.5$; $k_3 = f(0.5,\ 1 + 0.5\cdot1.5) = 1.75$; $k_4 = f(1,\ 1 + 1\cdot1.75) = 2.75$. Then $y_1 = 1 + \tfrac{1}{6}(1 + 3 + 3.5 + 2.75) = 1 + \tfrac{10.25}{6} = \tfrac{65}{24} = 2.708333$. In a *single* step of $h = 1$, RK4 lands within $0.01$ of the true $e$ — an error seventy times smaller than Euler's, for four evaluations of $f$ instead of one.

Here is why it is so good. For $y' = y$, RK4's one-step formula works out to $y_1 = y_0\left(1 + h + \tfrac{h^2}{2} + \tfrac{h^3}{6} + \tfrac{h^4}{24}\right)$ — the Taylor series of the exact multiplier $e^h$, correct through the $h^4$ term. RK4 reproduces the true solution's expansion up to fourth order and only errs at fifth, which is precisely what "fourth-order global accuracy" means. The whole ladder of methods, all on $y' = y$ with $h = 1$, tells the story at a glance:

   Euler      (1 eval /step)  ->  2.000000
   midpoint   (2 evals/step)  ->  2.500000
   RK4        (4 evals/step)  ->  2.708333
   exact  e                   ->  2.718282

⚡ Performance Note: RK4 costs four RHS evaluations per step versus Euler's one, so a naive glance says it is "4× slower per step." That is the wrong comparison. Because RK4 is fourth order, reaching a target accuracy $\varepsilon$ needs about $\varepsilon^{-1/4}$ steps, while Euler needs about $\varepsilon^{-1}$. For a modest accuracy of $10^{-6}$, Euler needs on the order of a million steps and RK4 on the order of thirty — so even at four evaluations each, RK4 does thousands of times less total work. Higher order is not a luxury; for smooth problems it is the cheapest accuracy you can buy.

📜 From History: The method is older than the computers it now runs on. Carl Runge published the core idea in 1895 and Wilhelm Kutta gave the fourth-order formula its modern form in 1901 — decades before there was any machine to run it on, as a way to integrate orbits and trajectories by hand. When the first electronic computers arrived in the 1950s, RK4 was waiting for them, and the earliest Fortran scientific libraries encoded it almost verbatim. The four lines you just typed are, coefficient for coefficient, what Kutta wrote down 125 years ago.

🔄 Check Your Understanding. 1. Write out the four RK4 stages from memory. What are the final weights, and why do they sum to one? 2. RK4 uses four evaluations of $f$ per step; Euler uses one. Why is RK4 nonetheless far cheaper to reach six-digit accuracy? 3. You suspect a colleague's "RK4" routine has a typo. What single numerical experiment exposes it?

Answers

1. $k_1 = f(t,y)$; $k_2 = f(t{+}\tfrac h2, y{+}\tfrac h2 k_1)$; $k_3 = f(t{+}\tfrac h2, y{+}\tfrac h2 k_2)$; $k_4 = f(t{+}h, y{+}h k_3)$; weights $\tfrac16,\tfrac13,\tfrac13,\tfrac16$. They sum to one so that a constant slope is integrated exactly (the step advances by $h\cdot\text{slope}$). 2. Fourth order means error $\sim h^4$, so a given accuracy needs dramatically fewer, larger steps — thousands of times fewer than first-order Euler despite 4× the per-step cost. 3. Halve $h$ and check the error ratio: true RK4 gives about $16\times$ ($2^4$); a mistyped stage collapses it to $\approx 4$ (second order) or $\approx 8$ (third).

23.3 Adaptive Step-Size Control

So far the step size $h$ has been a constant we picked in advance. That is wasteful. Real solutions have calm stretches, where a large step would be perfectly accurate, and violent stretches — a close orbital approach, a sudden reaction — where the same step would be wildly wrong. A fixed $h$ must be small enough for the worst moment, so it squanders effort everywhere else. The cure is to let the integrator measure its own error and choose its own step.

Definition (adaptive step size). An adaptive integrator estimates the local error of each tentative step and compares it to a user-specified tolerance. If the estimate exceeds the tolerance the step is rejected and retried with a smaller $h$; if it is comfortably under, the step is accepted and the next $h$ is grown. The step size tracks the difficulty of the solution, spending small steps only where they are needed.

The one new ingredient is an estimate of the local error, since we do not know the true solution. Two standard ways to get one:

Step doubling (Richardson). Take the step two ways — once with size $h$, and once as two steps of size $h/2$ — and compare. The two-half-step result is more accurate, and for a method of order $p$ the difference between them estimates the error of the finer result as $$\text{err} \approx \frac{y_{\text{half}} - y_{\text{full}}}{2^{p} - 1}.$$ This is the same Richardson-extrapolation idea you saw applied to derivatives and integrals in Chapter 22: combine two approximations at different resolutions to expose the leading error term.

Embedded pairs. More efficient, and what production code uses: an embedded Runge-Kutta pair computes two solutions of different orders from the same set of stages, so the error estimate is nearly free. The Runge-Kutta-Fehlberg method (RKF45) and the Dormand-Prince pair (the dopri5 behind MATLAB's ode45 and SciPy's default solve_ivp) are the famous examples: five or six stages yield both a fourth- and a fifth-order estimate, and their difference is the local error.

Once you have an error estimate err and a tolerance tol, the controller sets the next step from the order-$p$ scaling of the error, with a safety factor $S \approx 0.9$ to stay conservative: $$h_{\text{new}} = h \cdot S \cdot \left(\frac{\text{tol}}{\text{err}}\right)^{1/(p+1)}.$$ If err > tol this shrinks $h$ and the step is retried; if err < tol it grows $h$ for the next step.

The mechanics are clearest in the smallest honest example: one trial step of an adaptive Euler integrator ($p = 1$) on our test problem, using step doubling to estimate the error and the controller to propose the next step.

! adaptive_step.f90 -- one adaptive trial step (Euler + step doubling) for y' = y
program adaptive_trial
  implicit none
  integer,  parameter :: dp = selected_real_kind(15, 307)
  real(dp), parameter :: safety = 0.9_dp, tol = 0.01_dp
  integer,  parameter :: p = 1                     ! Euler is first order
  real(dp) :: t, y, h, y_full, y_half, err, h_new

  t = 0.0_dp;  y = 1.0_dp;  h = 0.5_dp

  y_full = y + h * f(t, y)                          ! one full step of size h
  y_half = y + 0.5_dp*h * f(t, y)                   ! two steps of size h/2
  y_half = y_half + 0.5_dp*h * f(t + 0.5_dp*h, y_half)

  err   = abs(y_half - y_full) / real(2**p - 1, dp) ! Richardson error estimate
  h_new = h * safety * (tol / err) ** (1.0_dp/real(p + 1, dp))

  print '(a, f10.6)', 'full step  y_full = ', y_full
  print '(a, f10.6)', 'two halves y_half = ', y_half
  print '(a, f10.6)', 'error estimate    = ', err
  print '(a, f10.6)', 'proposed next h   = ', h_new
contains
  pure function f(t, y) result(dydt)
    real(dp), intent(in) :: t, y
    real(dp)             :: dydt
    dydt = y
  end function f
end program adaptive_trial

$ gfortran -std=f2018 -Wall adaptive_step.f90 -o adapt && ./adapt
full step  y_full =   1.500000
two halves y_half =   1.562500
error estimate    =   0.062500
proposed next h   =   0.180000

By hand: the full step gives $y_{\text{full}} = 1 + 0.5(1) = 1.5$; the two half-steps give $1 + 0.25(1) = 1.25$ then $1.25 + 0.25(1.25) = 1.5625$. The error estimate is $(1.5625 - 1.5)/(2^1 - 1) = 0.0625$, which exceeds our tolerance of $0.01$, so this step would be rejected. The controller proposes $h_{\text{new}} = 0.5 \cdot 0.9 \cdot (0.01/0.0625)^{1/2} = 0.45 \cdot 0.4 = 0.18$, and we retry from $t = 0$ with the smaller step. (As a bonus, the Richardson-extrapolated value $y_{\text{half}} + 0.0625 = 1.625$ is far closer to the true $e^{0.5} = 1.6487$ than either raw step — free accuracy from the same two computations.)

🧩 Try It Yourself: Before reading on, predict what happens if you set tol = 0.1 instead of 0.01. Does this step get accepted or rejected, and does the next $h$ grow or shrink? Work it out with the formula, then change the one line and compile. (Answer: $\text{tol}/\text{err} = 0.1/0.0625 = 1.6$, so the step is accepted and $h$ grows to $0.5 \cdot 0.9 \cdot \sqrt{1.6} \approx 0.569$.)

🔗 Connection: You do not have to write adaptive Runge-Kutta from scratch for production work — the Fortran ecosystem provides it. ODEPACK's lsoda (a Tier-1 classic still in heavy use), the fortran-lang stdlib's ODE routines, and wrappers around Sundials all implement adaptive, order-aware stepping with decades of hardening. Write your own once to understand it, as we just did; reach for a battle-tested library when correctness and robustness matter. This is the ecosystem argument of Chapter 16 applied to ODEs.

23.4 Systems of ODEs and the Method of Lines

Real problems are never a single scalar equation. A planet moving in a plane has a position and a velocity, each with two components — four numbers evolving together. A chemical reactor has a concentration for every species. A predator-prey model couples two populations. The state is a vector, and the RHS returns a vector of rates: $$\frac{d\mathbf{y}}{dt} = \mathbf{f}(t, \mathbf{y}), \qquad \mathbf{y} \in \mathbb{R}^{m}.$$

Definition (first-order system). A first-order system of ODEs evolves a vector state $\mathbf{y} = (y_1, \dots, y_m)$ under a vector RHS $\mathbf{f}$. Crucially, any higher-order ODE can be rewritten as a first-order system by introducing the derivatives as new state variables. The second-order equation $y'' = g(t, y, y')$ becomes the pair $y_1' = y_2,\ y_2' = g(t, y_1, y_2)$ with $y_1 = y$ and $y_2 = y'$. This reduction is universal: master first-order systems and you can integrate any ODE, of any order.

The reduction is the key skill. Take the simple harmonic oscillator $y'' + y = 0$ — a mass on a spring, a pendulum at small amplitude. Set $y_1 = y$ (position) and $y_2 = y'$ (velocity); then $$y_1' = y_2, \qquad y_2' = -y_1,$$ with exact solution $y_1(t) = \cos t,\ y_2(t) = -\sin t$ for the start $(1, 0)$.

Here Fortran's arrays earn their keep. RK4 for a system is the same four stages, but each $k$ is now a vector and every operation is a whole-array operation from Chapter 5 — no inner loops, the code reads like the math:

! example-03-system.f90 -- one RK4 step of the harmonic oscillator y'' + y = 0
module sys_mod
  implicit none
  integer, parameter :: dp = selected_real_kind(15, 307)

  abstract interface
    function rhs_sys(t, y) result(dydt)
      import :: dp
      real(dp), intent(in) :: t
      real(dp), intent(in) :: y(:)
      real(dp)             :: dydt(size(y))
    end function rhs_sys
  end interface

contains

  function rk4_sys(f, t, y, h) result(y_next)
    procedure(rhs_sys)   :: f
    real(dp), intent(in) :: t, h, y(:)
    real(dp)             :: y_next(size(y))
    real(dp)             :: k1(size(y)), k2(size(y)), k3(size(y)), k4(size(y))
    k1 = f(t,            y)
    k2 = f(t + 0.5_dp*h, y + 0.5_dp*h*k1)
    k3 = f(t + 0.5_dp*h, y + 0.5_dp*h*k2)
    k4 = f(t + h,        y + h*k3)
    y_next = y + (h/6.0_dp) * (k1 + 2.0_dp*k2 + 2.0_dp*k3 + k4)
  end function rk4_sys

  function oscillator(t, y) result(dydt)   ! y1' = y2, y2' = -y1
    real(dp), intent(in) :: t
    real(dp), intent(in) :: y(:)
    real(dp)             :: dydt(size(y))
    dydt = [ y(2), -y(1) ]
  end function oscillator
end module sys_mod

program run_system
  use sys_mod
  implicit none
  real(dp) :: y(2), yn(2), energy

  y  = [1.0_dp, 0.0_dp]                 ! position 1, velocity 0
  yn = rk4_sys(oscillator, 0.0_dp, y, 0.5_dp)
  energy = 0.5_dp * (yn(1)**2 + yn(2)**2)
  print '(a, 2f12.6)', 'RK4 state (u1, u2) = ', yn
  print '(a, 2f12.6)', 'exact (cos, -sin)  = ', cos(0.5_dp), -sin(0.5_dp)
  print '(a,  f12.6)', 'energy after step  = ', energy
  print '(a,  f12.6)', 'exact energy       = ', 0.5_dp
end program run_system

$ gfortran -std=f2018 -Wall example-03-system.f90 -o system && ./system
RK4 state (u1, u2) =     0.877604   -0.479167
exact (cos, -sin)  =     0.877583   -0.479426
energy after step  =     0.499895
exact energy       =     0.500000

One RK4 step of $h = 0.5$ tracks the true $(\cos 0.5, -\sin 0.5)$ to about four decimals — and the energy $\tfrac{1}{2}(y_1^2 + y_2^2)$, which the true oscillator conserves exactly at $0.5$, has drifted only to $0.499895$. That the code for a two-component system is barely longer than the scalar version is entirely due to Fortran's array semantics: k1 = f(t, y) computes a whole vector, y + 0.5*h*k1 adds whole vectors, and the compiler is free to vectorize all of it. This is the payoff we promised back in Chapter 5 — arrays are Fortran's superpower, and nowhere more visibly than here.

💡 Intuition — energy drift. RK4 is superbly accurate over any fixed interval, but notice the energy crept below $0.5$. Over millions of steps that slow bleed matters: an RK4 planet spirals gently into its star. For long-term orbital work, scientists switch to symplectic integrators (leapfrog, velocity-Verlet) that are less accurate step-to-step but conserve a nearby energy exactly, so the orbit stays an orbit forever. Accuracy over a step and fidelity over an eon are different goals — a theme we return to in §23.6.

The method of lines: a PDE is a giant system of ODEs

Now the idea that reframes the whole book. A partial differential equation involves derivatives in several variables — for the heat equation, in time and space. But suppose we discretize only the space, laying down a grid and replacing the spatial derivatives with the finite differences of Chapter 22, while leaving time continuous. Then the temperature at each grid point becomes its own function of time, and the PDE turns into one ODE per grid point — a coupled first-order system exactly like the ones above, only enormous.

Definition (method of lines). The method of lines solves a time-dependent PDE by discretizing all spatial derivatives on a grid while keeping time continuous, producing a large system of ODEs $\frac{d\mathbf{u}}{dt} = \mathbf{F}(\mathbf{u})$ — one equation for the value at each grid point — which is then handed to any ODE integrator from this chapter. The spatial discretization and the time integration become two independent choices.

Make it concrete with the one-dimensional heat equation $\dfrac{\partial u}{\partial t} = \alpha\dfrac{\partial^2 u}{\partial x^2}$. Put down grid points $x_i$ spaced $\Delta x$ apart and approximate the second derivative with the standard central difference from Chapter 22, $\dfrac{\partial^2 u}{\partial x^2}\Big|_i \approx \dfrac{u_{i-1} - 2u_i + u_{i+1}}{\Delta x^2}$. The PDE becomes, at each interior point, $$\frac{du_i}{dt} = \frac{\alpha}{\Delta x^{2}}\left(u_{i-1} - 2u_i + u_{i+1}\right).$$ That is a linear system of ODEs, $\dfrac{d\mathbf{u}}{dt} = A\mathbf{u}$, where $A$ is the tridiagonal second-difference matrix. Nothing new is required to solve it — it is a first-order system, and you already have RK4 for those.

🚪 Threshold Concept — a PDE is a giant system of ODEs. This is the idea that unifies Part V. Once you discretize space, the distinction between "ODE solver" and "PDE solver" largely dissolves: a PDE is a system of ODEs with one equation per grid point, and marching it forward in time is exactly the initial-value problem of this chapter, scaled up from two components to a million. The method you choose for time-stepping a PDE — Euler, RK4, implicit — is an ODE method, chosen with the ODE reasoning you just learned. Every finite-difference PDE code in existence, including the one you are building, is an ODE integrator wearing a spatial-discretization hat.

And here is where it lands on your own project. Apply the simplest ODE method — Euler — to that system, and one step reads $$\mathbf{u}^{n+1} = \mathbf{u}^{n} + \Delta t\,\frac{\alpha}{\Delta x^2}\left(\text{second difference of }\mathbf{u}^{n}\right).$$ That is precisely the explicit heat-stepping scheme you will formalize in Chapter 24. The heat solver's time loop, which has felt like its own special construction since Chapter 4, turns out to be Euler's method in disguise — and that means RK4 could replace it verbatim for a higher-order-in-time solver. We cash this out in the Project Checkpoint below.

🔗 Connection — where the CFL condition comes from. The eigenvalues of that tridiagonal matrix $A$ are all real and negative, the most negative near $-4\alpha/\Delta x^2$. As you will see in §23.5, an explicit method is stable only when $\Delta t$ times the most negative eigenvalue stays inside a bounded region — for Euler, when $\Delta t\,(4\alpha/\Delta x^2) \le 2$, i.e. $\Delta t \le \Delta x^2/(2\alpha)$. That inequality is the CFL condition of Chapter 24, and now you know it is not a PDE mystery: it is the absolute-stability limit of an ODE method applied to the method-of-lines system. Shrink $\Delta x$ to resolve space and $\Delta t$ must shrink like $\Delta x^2$ to stay stable — the tyranny that motivates implicit stepping.

🔄 Check Your Understanding. 1. Rewrite the second-order ODE $y'' + 3y' + 2y = 0$ as a first-order system. 2. In the method of lines, which derivatives are discretized and which are left continuous? 3. What familiar heat-solver scheme do you get by applying Euler's method to the method-of-lines system?

Answers

1. Let $y_1 = y,\ y_2 = y'$; then $y_1' = y_2$ and $y_2' = -3y_2 - 2y_1$. 2. Spatial derivatives are discretized on the grid (finite differences); the time derivative is left continuous, giving one ODE per grid point. 3. The explicit forward-time, centered-space heat scheme — $\mathbf{u}^{n+1} = \mathbf{u}^n + \Delta t\,\alpha\,(\text{second difference})$ — the default explicit stepping of Chapter 24.

23.5 Stiffness and Implicit Methods — a Preview

There is a class of problems where everything above quietly fails, and it is common enough that you must recognize it: stiff problems. A stiff system has processes on wildly different time scales — a fast transient that dies almost instantly alongside a slow evolution you actually care about — and explicit methods like Euler and RK4 are forced to take absurdly small steps not for accuracy but for stability.

Definition (stiffness). A system of ODEs is stiff when it contains dynamics on very different time scales, so that an explicit integrator must use a step size far smaller than the accuracy of the solution would require — limited instead by stability. Symptomatically: the solution looks smooth and boring, yet your explicit solver either blows up or crawls. Stiffness is the rule, not the exception, in chemical kinetics, combustion, electronics, and the finely resolved PDEs of §23.4.

The cleanest way to see it is the linear test equation $y' = \lambda y$ with $\lambda < 0$, whose true solution $y = y_0 e^{\lambda t}$ decays smoothly to zero. Apply explicit Euler: $y_{n+1} = (1 + h\lambda)\,y_n$. The numerical solution decays only if the amplification factor satisfies $|1 + h\lambda| \le 1$, which for real $\lambda < 0$ means $h \le 2/|\lambda|$. Take $\lambda = -1$ and watch the step size decide the fate of a solution that should just fade away:

   y' = -y,  y(0) = 1        one explicit Euler step, y1 = 1 + h*(-1)
   h = 1 :  y1 =  0          (decays to the fixed point -- fine)
   h = 2 :  y1 = -1          (marginal: oscillates +/-1 forever)
   h = 3 :  y1 = -2          (|amplification| > 1 -- blows up: -2, 4, -8, ...)

For $\lambda = -1$ the limit is $h \le 2$, hardly a burden. But a stiff system has some $\lambda \approx -10^{6}$, and then explicit Euler demands $h \le 2\times10^{-6}$ for the entire run, even during the eons when nothing is happening — millions of tiny steps to cross a smooth stretch. RK4 helps only marginally: its real-axis stability limit is about $h|\lambda| \le 2.8$ instead of $2$, a 40% reprieve, not a solution. Explicit methods are simply the wrong tool for stiff problems.

Definition (implicit method). An implicit method defines the new state in terms of the RHS evaluated at the new state itself. Backward (implicit) Euler is $$y_{n+1} = y_n + h\,f(t_{n+1},\, y_{n+1}),$$ with $y_{n+1}$ appearing on both sides. You cannot simply evaluate it; you must solve an equation for $y_{n+1}$ each step. In exchange you get vastly better stability.

Apply backward Euler to $y' = \lambda y$: $y_{n+1} = y_n + h\lambda\,y_{n+1}$, so $y_{n+1} = y_n/(1 - h\lambda)$. For $\lambda < 0$ the denominator $1 - h\lambda$ exceeds one for every $h > 0$, so $|y_{n+1}| < |y_n|$ always — the numerical solution decays for any step size whatsoever. This property is called A-stability, and it is why implicit methods own the stiff regime. On the same $y' = -y$, backward Euler with $h = 3$ gives $y_1 = 1/(1 - 3(-1)) = 1/4 = 0.25$: still decaying, still sane, where explicit Euler had already exploded.

The price is the equation solve. For a linear system, $y_{n+1} = y_n/(1 - h\lambda)$ generalizes to solving a linear system $(I - h A)\mathbf{y}_{n+1} = \mathbf{y}_n$ every step — a job for the LAPACK solvers of Chapter 21, which is exactly why that chapter's optional Project Checkpoint set up a backward-Euler heat step as a tridiagonal solve. For a nonlinear RHS you need Newton's method inside each step. Either way, an implicit step costs far more than an explicit one — but for a stiff problem it lets you take steps thousands of times larger, and wins by a landslide. This is a preview; the explicit-versus-implicit choice for the heat equation is developed in full in Chapter 24.

🐛 Find the Bug. A colleague integrates a stiff chemical network with RK4, halves the step until it stops blowing up, and reports that the code "works but is unbearably slow — hours for one second of simulated time." Their RHS and coefficients are all correct. What is actually wrong, and what should they change?

Answer

Nothing is wrong with the code — the method is wrong for the problem. The tiny step is forced by stability, not accuracy: a fast-decaying mode with a large negative $\lambda$ is pinning $h$ below $\approx 2.8/|\lambda|$ even though that mode died in the first microsecond. RK4 is explicit and cannot escape this. The fix is an implicit or specialized stiff solver (backward Euler, BDF, or a library like ODEPACK's lsoda, which detects stiffness and switches methods automatically), which is A-stable and can take large steps through the smooth regime.

23.6 Applications: Orbits, Reactions, and Populations

The methods of this chapter are general, but they come alive in specific problems. Three families cover an enormous share of scientific ODE work, and each teaches something the others do not.

Orbital motion. Newton's law of gravitation for a body orbiting a fixed mass gives $\ddot{\mathbf{r}} = -GM\,\mathbf{r}/|\mathbf{r}|^3$. In a plane this is a second-order vector equation; reduce it in the usual way to a first-order system of four — $x,\ y,\ v_x,\ v_y$ — and integrate. The lesson here is conservation: the true orbit conserves energy and angular momentum exactly, so those quantities are your error diagnostic. Watch the energy, and you will see RK4's slow drift (§23.4) and understand why long-duration ephemeris and molecular-dynamics codes reach for symplectic integrators instead. Integrating an orbit is the classic first project because the answer is a shape you can see — a closed ellipse if you did it right, a slow spiral if your method leaks energy.

Chemical kinetics. A reaction network is a system of ODEs in the species concentrations, with the RHS built from the reaction rates. A simple linear decay chain $A \to B \to C$ is $$[A]' = -k_1[A], \qquad [B]' = k_1[A] - k_2[B], \qquad [C]' = k_2[B],$$ which conserves total mass ($[A]+[B]+[C]$ is constant) — another built-in sanity check. Kinetics is also where you meet stiffness in the wild: real mechanisms have rate constants spanning many orders of magnitude (a radical that reacts in nanoseconds coupled to a product that forms over minutes), so production combustion and atmospheric-chemistry codes are built on the implicit stiff solvers of §23.5. The famous Robertson problem — three species, rate constants of $0.04$, $10^4$, and $3\times10^7$ — is the standard stiff-solver benchmark for exactly this reason.

Population dynamics. The logistic equation $N' = rN(1 - N/K)$ models growth that saturates at a carrying capacity $K$; it is a single nonlinear ODE with a known solution, ideal for verifying a solver. Couple two species and you get Lotka-Volterra predator-prey, $$x' = \alpha x - \beta x y, \qquad y' = \delta x y - \gamma y,$$ whose solutions are closed loops in the $(x, y)$ plane — populations oscillate forever, predator lagging prey. Like the orbit, it conserves a quantity, so integrating it well means keeping the loop closed; integrate it with plain Euler and the spurious energy gain spirals the populations outward, a vivid visual failure of a first-order method.

🔗 Connection: All three families are systems, all reduce to the first-order vector form of §23.4, and all are integrated by the very same rk4_sys routine you wrote — only the RHS function changes. That is the power of passing the RHS as a procedure argument (Chapter 6): one hardened integrator, any number of physical models. Swapping oscillator for a kepler or lotka_volterra function is the whole edit.

Project Checkpoint

This checkpoint reframes the heat solver, and it is the conceptual hinge of the whole project. You have been building a time loop since Chapter 4 without a name for what it is. Here is the name: your heat solver is an ODE integrator, and its time-stepping is Euler's method applied to the method-of-lines system.

Recall the setup from §23.4. Discretize the plate in space and each grid temperature becomes a function of time obeying $\dfrac{du_i}{dt} = \dfrac{\alpha}{\Delta x^2}(\text{second difference of } u)$ — a large first-order system $\mathbf{u}' = \mathbf{F}(\mathbf{u})$. Marching it with Euler is the explicit heat scheme. To make this unmistakable, here is a one-dimensional rod with three interior points, integrated by explicit Euler in time — with the heat RHS written as its own function, ready to be handed to any integrator in this chapter:

! project-checkpoint.f90 -- heat time-stepping IS Euler on the method-of-lines system
program mol_heat
  implicit none
  integer,  parameter :: dp = selected_real_kind(15, 307)
  integer,  parameter :: n = 3                 ! interior grid points
  real(dp), parameter :: alpha = 1.0_dp, dx = 1.0_dp, dt = 0.1_dp
  real(dp), parameter :: u_left = 1.0_dp, u_right = 0.0_dp  ! hot / cold ends
  real(dp) :: u(n)
  integer  :: step

  u = 0.0_dp                                   ! interior starts cold
  print '(a)', '1D heat by the method of lines (explicit Euler in time):'
  print '(a, 3f11.6)', 'step 0: ', u
  do step = 1, 2
    u = u + dt * heat_rhs(u)                   ! <-- one Euler step of the ODE system
    print '(a, i0, a, 3f11.6)', 'step ', step, ': ', u
  end do
contains
  function heat_rhs(u) result(dudt)            ! F(u): (alpha/dx^2) * second difference
    real(dp), intent(in) :: u(:)
    real(dp)             :: dudt(size(u)), c
    integer              :: i
    c = alpha / dx**2
    dudt(1) = c * (u_left  - 2.0_dp*u(1) + u(2))
    do i = 2, size(u) - 1
      dudt(i) = c * (u(i-1) - 2.0_dp*u(i) + u(i+1))
    end do
    dudt(size(u)) = c * (u(size(u)-1) - 2.0_dp*u(size(u)) + u_right)
  end function heat_rhs
end program mol_heat

$ gfortran -std=f2018 -Wall project-checkpoint.f90 -o mol_heat && ./mol_heat
1D heat by the method of lines (explicit Euler in time):
step 0:    0.000000   0.000000   0.000000
step 1:    0.100000   0.000000   0.000000
step 2:    0.180000   0.010000   0.000000

Trace it: with $\alpha = \Delta x = 1$, the rate at point 1 in step 1 is $u_{\text{left}} - 2u_1 + u_2 = 1 - 0 + 0 = 1$, so $u_1 \leftarrow 0 + 0.1(1) = 0.1$; points 2 and 3 see no gradient yet and stay $0$. In step 2 the heat has begun to spread: $u_1 \leftarrow 0.1 + 0.1(1 - 0.2) = 0.18$ and $u_2 \leftarrow 0 + 0.1(0.1) = 0.01$. Heat marches inward from the hot end, one point per step — diffusion, emerging from nothing but u = u + dt*heat_rhs(u).

Two things to carry into Chapter 24. First, the line u = u + dt * heat_rhs(u) is literally euler_step from §23.1 with a vector state — so swapping in rk4_sys from §23.4 gives a higher-order-in-time heat solver for free, the same spatial stencil integrated more accurately in time. Second, our $\Delta t = 0.1$ was not arbitrary: the stability limit here is $\Delta t \le \Delta x^2/(2\alpha) = 0.5$, and $0.1$ sits safely under it. Push $\Delta t$ past $0.5$ and this explicit scheme blows up — the CFL condition of §23.4, waiting for you in Chapter 24. Save this file as heat-solver/heat_mol.f90; it is the explicit core, and Chapter 24 makes it real in two dimensions.

Summary

This chapter turned the abstract "solve a differential equation" into concrete, compilable time-stepping, and connected it to the heat solver.

The two things to memorize: the four RK4 stages exactly (a mistyped coefficient silently destroys the order), and the one-line reframing that u = u + dt*rhs(u) is Euler's method — so the heat solver is an ODE integrator, and a better time integrator drops straight in.

Terms first defined here: initial-value problem, right-hand-side function, Euler's method, step size, local truncation error / global order (for integrators), Runge-Kutta, RK4, Butcher tableau, adaptive step size, first-order system, method of lines, stiffness, implicit method / backward Euler, A-stability.

Compile note: every example in this chapter builds with the plain gfortran -std=f2018 -Wall file.f90 -o exe — no external libraries. The RHS is passed as a procedure argument through an abstract interface, so one integrator serves any equation.

Spaced Review

Retrieval practice on Chapter 6 (procedures) and Chapter 22 (accuracy and finite differences) — the two chapters this one leans on hardest.

1. Our integrators take the RHS f as a dummy procedure declared procedure(rhs) :: f, and every real argument carries an intent. Why is intent(in) the right choice for the RHS's $t$ and $y$ arguments, and what does declaring it buy you? (Chapter 6)
2. The system RHS uses an assumed-shape array argument, real(dp), intent(in) :: y(:). What does the (:) mean, and why is it preferable to passing the size as a separate integer argument? (Chapter 6)
3. Marking the little test RHS pure (as in adaptive_step.f90) is a promise to the compiler. What exactly does pure forbid, and why does it help here? (Chapter 6)
4. RK4 is fourth-order and Euler is first-order. Using the Chapter 22 notion of order of accuracy, predict the error ratio when you halve $h$ for each method, and say how you would confirm the order numerically. (Chapter 22)
5. The method of lines replaces $\partial^2 u/\partial x^2$ with the central second difference $(u_{i-1} - 2u_i + u_{i+1})/\Delta x^2$. What is its order of accuracy in $\Delta x$, and how did Chapter 22 have you verify a finite-difference order in practice? (Chapter 22)

What's Next

You now hold every idea the heat solver needs, and one of them — the method of lines — has quietly turned the remaining PDE work into applied ODE integration. Chapter 24 cashes that in. It builds the finite-difference stencil for the two-dimensional heat equation in full, makes the explicit time-stepping you just previewed into the solver's real core, derives the CFL stability condition you glimpsed here as an ODE stability limit, and handles boundary conditions properly. The u = u + dt*rhs(u) line from this chapter's checkpoint becomes the beating heart of your simulation. Let's discretize the plate.