Systems of ODEs — RK4 for systems: Visual walkthrough

We assume you can add, multiply, and read a graph. That is all. Every new symbol is drawn before it is used.

Step 1 — What "solving an ODE" even means (the slope field)

WHAT. We have an unknown quantity that changes over time. Call the time $t$ and the quantity $y$ . An ordinary differential equation (ODE) is a rule that tells you the slope of $y$ at any point: $y^{'} = f (t, y) .$

Let us name every symbol on this line, right where it sits:

$y^{'}$ — the slope, i.e. how fast $y$ is rising or falling. Read it "y-prime".
$f (t, y)$ — a known formula that, given where you are ( $t$ ) and how high you are ( $y$ ), hands you back that slope.

WHY. We do not know $y$ itself yet — that is the whole puzzle. But we are handed its slope everywhere. So solving means: start at a known point and walk, always obeying the slope the formula gives.

PICTURE. Figure s01 draws the slope field for $y^{'} = y$ : at every point a tiny arrow shows which way to walk. The true solution (the smooth curve $e^{t}$ ) is the path that stays tangent to every arrow it passes.

Figure — Systems of ODEs — RK4 for systems

Step 2 — One step is really an area (the exact integral)

WHAT. Take one small time-step of width $h$ , from $t_{n}$ to $t_{n} + h$ . The exact change in $y$ over that step is $y (t_{n} + h) = y (t_{n}) + shaded area \int_{t_{n}}^{t_{n} + h} f (s, y (s)) d s .$

Term by term:

$h$ — the step width, a small chosen number like $0.1$ .
$\int_{t_{n}}^{t_{n} + h} (\dots) d s$ — the "add-up-all-the-slopes" symbol; it sums the slope $f$ across the whole step. The letter $s$ is just a stand-in for time as it sweeps from $t_{n}$ to $t_{n} + h$ .
The result is the area under the slope curve over that step.

WHY. "Rise = area under the slope curve" is the fundamental fact of calculus. If we could compute that area exactly, we would land exactly on the true curve. We cannot — because $y (s)$ inside $f$ is the very thing we don't know. So the entire game becomes: estimate this one area cleverly.

PICTURE. Figure s02 shades the area under the slope curve across one step. Landing accuracy = how well we guess this shaded region.

Step 3 — Simpson's idea: sample the slope at start, middle, end

WHAT. Simpson's Rule estimates an area by sampling the curve at three places — the two ends and the midpoint — and weighting the middle four times as heavily: $\int_{t_{n}}^{t_{n} + h} f d s \approx \frac{h}{6} (\times 1 f_{start} + \times 4 4 f_{mid} + \times 1 f_{end}) .$

Reading the weights: $1$ at the start, $4$ in the middle, $1$ at the end — total $6$ , which is exactly the $6$ dividing the outside.

WHY THIS TOOL AND NOT ANOTHER. Why Simpson, not just "slope at the start times $h$ " (that would be Euler)? Because a single starting slope ignores how the slope bends during the step. Simpson fits a parabola through three samples — it captures curvature, so it is dramatically more accurate for the same step size. That extra accuracy is what will earn RK4 its order-4 rating.

PICTURE. Figure s03 overlays a parabola on three sampled points and shows the $1, 4, 1$ weighting as bar heights.

Step 4 — Building the four slope estimates $k_{1}, k_{2}, k_{3}, k_{4}$

WHAT. We compute four slope samples, each using a better guess of where we are: $k_{1} = f (t_{n}, y_{n})$ $k_{2} = f (t_{n} + \frac{h}{2}, y_{n} + \frac{h}{2} k_{1})$ $k_{3} = f (t_{n} + \frac{h}{2}, y_{n} + \frac{h}{2} k_{2})$ $k_{4} = f (t_{n} + h, y_{n} + h k_{3})$

Symbol-by-symbol, what each guess is doing:

$k_{1}$ — slope at the start, using the state $y_{n}$ we already know.
$y_{n} + \frac{h}{2} k_{1}$ — an Euler half-step forward using $k_{1}$ : "if the start-slope held, here's the midpoint." $k_{2}$ = slope there.
$y_{n} + \frac{h}{2} k_{2}$ — a refined midpoint, now using the better slope $k_{2}$ . $k_{3}$ = slope there.
$y_{n} + h k_{3}$ — a full step using the good midpoint slope $k_{3}$ , landing near the end. $k_{4}$ = slope there.

WHY. Each estimate feeds the next: we bootstrap our way to trustworthy midpoint and endpoint slopes without ever knowing the true curve. $k_{2}$ and $k_{3}$ are two independent estimates of the same midpoint slope — averaging them cancels errors.

PICTURE. Figure s04 draws all four probe arrows on one step: $k_{1}$ at the left, the two midpoint probes stacked, $k_{4}$ at the right.

Step 5 — The weighted average that lands the step

WHAT. Combine the four slopes with Simpson's spirit — but Simpson had only one middle sample, and we have two ( $k_{2}, k_{3}$ ). So we split Simpson's weight $4$ evenly, giving each middle slope weight $2$ : $y_{n + 1} = y_{n} + \frac{h}{6} (\times 1 k_{1} + \times 2 2 k_{2} + \times 2 2 k_{3} + \times 1 k_{4}) .$

The weights $1 + 2 + 2 + 1 = 6$ match the divisor $6$ — this guarantees that if the slope were constant, we'd get exactly slope $\times h$ (no bias).

WHY. This precise averaging makes the error terms cancel out to fourth order — the leftover error per step is $O (h^{5})$ , tiny (halving $h$ shrinks it 32×). Accumulated over a whole interval the global error is $O (h^{4})$ ; see Local vs Global Truncation Error.

PICTURE. Figure s05 shows the four slope arrows fanned out and their weighted average arrow (in orange) pointing to the landing point.

Step 6 — Now let $y$ become a vector $y$

WHAT. Real problems track several quantities at once (position and velocity, prey and predator). Stack them into one list — a vector: $y = y_{1} y_{2} ⋮ y_{n}, f (t, y) = f_{1} (t, y_{1}, \dots, y_{n}) ⋮ f_{n} (t, y_{1}, \dots, y_{n}) .$

The bold $y$ means "the whole list at once." Crucially, each slope-formula $f_{i}$ may read every component — that is coupling: the components hold hands.

WHY THE DERIVATION IS FREE. Look back at Steps 4–5: we only ever added things and multiplied by numbers like $\frac{h}{2}$ . Those two operations work identically on a list of numbers as on a single number. So we bold every symbol and change nothing else: $k_{1} = f (t_{n}, y_{n}), k_{2} = f (t_{n} + \frac{h}{2}, y_{n} + \frac{h}{2} k_{1}), \dots$ $y_{n + 1} = y_{n} + \frac{h}{6} (k_{1} + 2 k_{2} + 2 k_{3} + k_{4})$ Each $k_{i}$ is now a vector of length $n$ .

PICTURE. Figure s06 shows two coupled trajectories $y_{1} (t), y_{2} (t)$ advancing in lockstep — one step of the vector recipe moving both at once.

Step 7 — The degenerate & edge cases (never leave a gap)

WHAT / WHY / PICTURE for the corners you must handle:

One equation ( $n = 1$ ). The vector collapses to a scalar; bold marks vanish; you recover Step 5 exactly. RK4-for-systems contains scalar RK4 as a special case.
Uncoupled system (each $f_{i}$ depends only on its own $y_{i}$ ). Then the equations don't hold hands and you could run separate RK4s — but the vector recipe still works and costs nothing extra. Coupling is the general case; independence is the lucky one.
The synchronization trap (the one real mistake). Because $f_{1}$ can read $y_{2}$ , you must finish the entire vector $k_{1}$ before forming any intermediate state. If you fully finish RK4 on $y_{1}$ first, then $y_{1}$ 's $k_{2}$ reads a stale, un-advanced $y_{2}$ — silent, wrong answers.

Figure s07 contrasts the RIGHT order (all components advance together at each stage) against the WRONG order (one component races ahead, feeding stale data to the other).

Worked check — reproduce the parent's Example 1

For $y_{1}^{'} = y_{2}, y_{2}^{'} = - y_{1}$ , $y (0) = (1, 0)$ , $h = 0.1$ , one step gives (parent's numbers): $y_{1, 1} \approx 0.995004, y_{2, 1} \approx - 0.099833,$ matching $cos (0.1) = 0.995004$ and $- sin (0.1) = - 0.099833$ . The VERIFY block below recomputes both from scratch.

The one-picture summary

Figure s08 compresses the whole walkthrough: the exact-area problem (Step 2), the three-sample Simpson estimate split into four probes (Steps 3–4), the $1, 2, 2, 1$ weighted landing (Step 5), and the vector bolding that makes it all work for coupled systems (Step 6).

Recall Feynman retelling (plain words)

We want to know where a moving thing will be one moment from now. The true answer is the area under its slope curve for that moment — but we can't see that area because it depends on where the thing goes, which is what we're solving for. So we cheat cleverly: we peek at the slope four times. Once at the very start. Twice in the middle — the second peek smarter than the first because it uses what the first peek learned. Once near the end using the best middle guess. Then we average the four peeks, counting the two middle ones double, and divide by six. That average slope, times the step, tells us how far to move. When several things move together holding hands, we do the exact same peeking — but every peek must move all of them at once, never one alone, or they'd be reading each other's out-of-date positions. That is RK4 for systems: four synchronized peeks, one weighted average, over and over.

Connections

Parent topic
RK4 for a single ODE — the scalar recipe we bolded.
Simpson's Rule — the area-sampling RK4 imitates.
Reducing higher-order ODEs to first-order systems
Euler's method for systems — the one-peek cousin.
Local vs Global Truncation Error — why $O (h^{5}) \to O (h^{4})$ .
Stiff systems and stability — where explicit RK4 struggles.