Worked examples — Systems of ODEs — RK4 for systems

Before any example, one reminder of the machinery we reuse every time.

Here "vector" just means an ordered list of numbers, one per unknown. If we track two quantities $y_1,y_2$, then $\mathbf{y}=(y_1,y_2)$ is a point in a plane, and $\mathbf{f}$ hands back another such pair — the direction and speed each quantity is currently moving. That is all the notation means.

---

The scenario matrix

Every problem this topic can throw at you falls into one of these cells. The worked examples below are each tagged with the cell(s) they cover, so that together they fill the whole grid.

Cell	What makes it different	Example
A. Coupled, mixed signs	$f_1,f_2$ each depend on the other variable, opposite signs (rotation)	Ex 1
B. Decoupled system	equations don't reference each other — a sanity check	Ex 2
C. Reduced 2nd-order	a single higher-order ODE turned into a system	Ex 3
D. Nonlinear word problem	products of unknowns, real units (predator–prey)	Ex 4
E. Degenerate / steady state	initial state where $\mathbf{f}=0$ (nothing should move)	Ex 5
F. All-positive growth	both variables grow, positive slopes throughout	Ex 6
G. Limiting / stiff twist	a component decays so fast that step size matters	Ex 7
H. Non-autonomous	$\mathbf{f}$ depends on $t$ explicitly, not just on $\mathbf{y}$	Ex 8

The columns of "sign behaviour" (both slopes negative, both positive, mixed, zero) and the rows of "structure" (coupled / decoupled / reduced / nonlinear / stiff / time-dependent) are all touched. The figure below plots four of these cells as motions in state space, so you can see how different the scenarios are before we grind through the arithmetic.

Read the figure like a map: the blue arc (Cell A, Ex 1) is a point curving on a circle — mixed signs make it rotate; the orange arrow (Cell B, Ex 2) shows one coordinate shrinking while the other grows — a decoupled system; the green curve (Cell F, Ex 6) is both coordinates racing outward together — pure growth; the red dot (Cell E, Ex 5) sits frozen at a fixed point — nothing moves. Keep this picture in mind as each example below traces one of these paths.

---

Ex 1 — Cell A: coupled rotation, mixed signs

Forecast: these are the equations of a point spinning on a circle (position and velocity feeding each other with opposite signs). After a small step $y_1$ should dip slightly below $1$ and $y_2$ should go slightly negative. Since RK4 is order 4, expect ~4 correct digits.

Here $\mathbf{f}(t,\mathbf{y})=(y_2,-y_1)$. This is the blue arc in the figure above.

1. ${\mathbf{k}}_1=(y_2,-y_1)=(0,-1)$. Why this step? The slope at the current state $(1,0)$ — this is where $\mathbf{y}$ is heading right now.
2. Midpoint state $=\mathbf{y}+\frac{h}{2}{\mathbf{k}}_1=(1,0)+0.1(0,-1)=(1,-0.1)$, so ${\mathbf{k}}_2=(-0.1,-1)$. Why this step? Re-measure the slope at the predicted middle of the step — both components use the advanced state.
3. $\mathbf{y}+\frac{h}{2}{\mathbf{k}}_2=(1,0)+0.1(-0.1,-1)=(0.99,-0.1)$, so ${\mathbf{k}}_3=(-0.1,-0.99)$. Why this step? A refined midpoint slope, correcting ${\mathbf{k}}_2$.
4. $\mathbf{y}+h{\mathbf{k}}_3=(1,0)+0.2(-0.1,-0.99)=(0.98,-0.198)$, so ${\mathbf{k}}_4=(-0.198,-0.98)$. Why this step? The end-of-interval slope.
5. Combine with weights $1,2,2,1$ over $6$: $y_{1,1}=1+\frac{0.2}{6}(0+2(-0.1)+2(-0.1)+(-0.198))=1+\frac{0.2}{6}(-0.598)\approx \mathrm{0.980067.}$ $y_{2,1}=0+\frac{0.2}{6}(-1+2(-1)+2(-0.99)+(-0.98))=\frac{0.2}{6}(-5.96)\approx -\mathrm{0.198667.}$ Why this step? We now blend the four slope vectors into one best estimate of the step's average velocity (middle slopes doubled — Simpson weighting) and move $\mathbf{y}$ by that average times $h$. This is the actual advance of the state.

Verify: $\cos 0.2=0.980067$ (matches to ~5 decimals), $-\sin 0.2=-0.198669$ (matches to ~4 decimals). Agreement in one step. Forecast confirmed.

---

Ex 2 — Cell B: decoupled sanity check

Forecast: the two equations never mention each other. RK4 on a decoupled system must reproduce, component-by-component, exactly what scalar RK4 gives on each equation separately. $y_1$ shrinks (negative slope), $y_2$ grows fast (slope $2y_2$). This is the orange arrow in the figure above.

Here $\mathbf{f}=(-y_1,2y_2)$.

1. ${\mathbf{k}}_1=(-1,2)$. Why? Slope at $(1,1)$.
2. State $=(1,1)+0.25(-1,2)=(0.75,1.5)$, ${\mathbf{k}}_2=(-0.75,3)$. Why? Midpoint slope.
3. State $=(1,1)+0.25(-0.75,3)=(0.8125,1.75)$, ${\mathbf{k}}_3=(-0.8125,3.5)$. Why? Refined midpoint.
4. State $=(1,1)+0.5(-0.8125,3.5)=(0.59375,2.75)$, ${\mathbf{k}}_4=(-0.59375,5.5)$. Why? End slope.
5. Combine: $y_{1,1}=1+\frac{0.5}{6}(-1+2(-0.75)+2(-0.8125)+(-0.59375))=1+\frac{0.5}{6}(-4.71875)\approx \mathrm{0.606771.}$ $y_{2,1}=1+\frac{0.5}{6}(2+2(3)+2(3.5)+5.5)=1+\frac{0.5}{6}(20.5)\approx \mathrm{2.708333.}$ Why this step? Each component gets its own $1,2,2,1$ blend — since the equations are decoupled, this is literally scalar RK4 done twice in parallel.

Verify: $e^{-0.5}=0.606531$ (error $\approx 2.4\times 10^{-4}$), $e^1=2.718282$ (error $\approx 0.01$). The larger growth rate ($2y_2$) makes a bigger step per unit $t$, so its error is larger — expected. Both are the same numbers you'd get running scalar RK4 on each line alone.

---

Ex 3 — Cell C: reduce a 2nd-order ODE

Forecast: this is the reduction trick. With $y(0)=y^{\prime }(0)=1$ the exact solution is $e^t$, so $u_1$ should track $e^{0.5}=1.6487$.

Set $u_1=y,u_2=y^{\prime }$. Then $u_1^{\prime }=u_2,u_2^{\prime }=u_1,$ so $\mathbf{f}=(u_2,u_1)$, starting at $\mathbf{u}=(1,1)$.

1. ${\mathbf{k}}_1=(1,1)$. Why? Slope at $(1,1)$: $u_1^{\prime }=u_2=1$, $u_2^{\prime }=u_1=1$.
2. State $=(1,1)+0.25(1,1)=(1.25,1.25)$, ${\mathbf{k}}_2=(1.25,1.25)$. Why? Midpoint slope — swap the components.
3. State $=(1,1)+0.25(1.25,1.25)=(1.3125,1.3125)$, ${\mathbf{k}}_3=(1.3125,1.3125)$. Why? Refined midpoint.
4. State $=(1,1)+0.5(1.3125,1.3125)=(1.65625,1.65625)$, ${\mathbf{k}}_4=(1.65625,1.65625)$. Why? End slope.
5. Combine (both components identical here since $u_1=u_2$ throughout): $u_{1,1}=1+\frac{0.5}{6}(1+2(1.25)+2(1.3125)+1.65625)=1+\frac{0.5}{6}(7.78125)\approx \mathrm{1.648438.}$ Why this step? Blend the four slopes with the fixed $1,2,2,1$ weights to advance $u_1$ (and $u_2$) by one step; the reduction changed nothing about this averaging.

Verify: $e^{0.5}=1.648721$; error $\approx 2.8\times 10^{-4}$. The 2nd-order ODE was solved by the same systems recipe — the reduction cost us nothing.

---

Ex 4 — Cell D: nonlinear predator–prey word problem

Forecast: the products $RF$ make this nonlinear — the coupling is multiplicative. At the start $R^{\prime }=2-2=0$ and $F^{\prime }=2-1=1$: rabbits momentarily steady, foxes rising. So expect $R\approx 2$ (nearly flat) and $F$ up to about $1.1$.

Here $\mathbf{f}=(R-RF,RF-F)$, state $(2,1)$.

1. ${\mathbf{k}}_1=(2-2\cdot 1,2\cdot 1-1)=(0,1)$. Why? Current slopes.
2. State $=(2,1)+0.05(0,1)=(2,1.05)$; ${\mathbf{k}}_2=(2-2\cdot 1.05,2\cdot 1.05-1.05)=(-0.1,1.05)$. Why? Midpoint slope, both populations advanced.
3. State $=(2,1)+0.05(-0.1,1.05)=(1.995,1.0525)$; $RF=1.995\cdot 1.0525=2.09974$, so ${\mathbf{k}}_3=(1.995-2.09974,2.09974-1.0525)=(-0.104738,1.047238)$. Why? Refined midpoint.
4. State $=(2,1)+0.1(-0.104738,1.047238)=(1.989526,1.104724)$; $RF=1.989526\cdot 1.104724=2.197931$, so ${\mathbf{k}}_4=(1.989526-2.197931,2.197931-1.104724)=(-0.208405,1.093207)$. Why? End slope.
5. Combine: $R_1=2+\frac{0.1}{6}(0+2(-0.1)+2(-0.104738)+(-0.208405))=2+\frac{0.1}{6}(-0.617881)\approx \mathrm{1.989702.}$ $F_1=1+\frac{0.1}{6}(1+2(1.05)+2(1.047238)+1.093207)=1+\frac{0.1}{6}(6.287683)\approx \mathrm{1.104795.}$ Why this step? Even though $\mathbf{f}$ is nonlinear, the averaging step is unchanged — we still weight the four slope vectors $1,2,2,1$ over $6$. Nonlinearity only affects how each ${\mathbf{k}}_i$ was computed, not how they are blended.

Verify: rabbits $\approx 1.9897$ thousand (barely dropped, as forecast — they were near their peak), foxes $\approx 1.1048$ thousand (rose, as forecast). Units are thousands of animals, consistent throughout; no negative populations appeared.

---

Ex 5 — Cell E: degenerate steady state

Forecast: at $(1,1)$ we get $R^{\prime }=1-1=0$ and $F^{\prime }=1-1=0$ — the field is exactly zero. A steady state means nothing should move: the answer must be $(1,1)$ again, to machine precision. This is the degenerate case that catches sloppy code — the red dot in the figure above.

1. ${\mathbf{k}}_1=(1-1,1-1)=(0,0)$. Why? Slope at the fixed point is zero by definition.
2. State $=(1,1)+0.15(0,0)=(1,1)$; ${\mathbf{k}}_2=(0,0)$. Why? We never left $(1,1)$, so the slope is still zero.
3. State $=(1,1)$; ${\mathbf{k}}_3=(0,0)$. Why? Same reasoning.
4. State $=(1,1)$; ${\mathbf{k}}_4=(0,0)$. Why? Same.
5. Combine: ${\mathbf{y}}_1=(1,1)+\frac{0.3}{6}(0+0+0+0)=(1,1)$. Why this step? The $1,2,2,1$ blend of four zero vectors is zero, so the average velocity is zero and the state does not advance — exactly what an equilibrium demands.

Verify: the state is unchanged — RK4 respects equilibria exactly (since every ${\mathbf{k}}_i=0$). If your program had "advanced one component alone," it might have nudged off the fixed point; the synchronized recipe does not.

---

Ex 6 — Cell F: all-positive growth

Forecast: with $\mathbf{y}(0)=(1,1)$ and the field pointing "each component grows at the other's value," both slopes are positive throughout — pure exponential growth. Expect both to reach $e^{0.4}=1.4918$. This is the green curve in the figure above.

Here $\mathbf{f}=(y_2,y_1)$, state $(1,1)$.

1. ${\mathbf{k}}_1=(1,1)$. Why? Slope at start.
2. State $=(1,1)+0.2(1,1)=(1.2,1.2)$, ${\mathbf{k}}_2=(1.2,1.2)$. Why? Midpoint.
3. State $=(1,1)+0.2(1.2,1.2)=(1.24,1.24)$, ${\mathbf{k}}_3=(1.24,1.24)$. Why? Refined midpoint.
4. State $=(1,1)+0.4(1.24,1.24)=(1.496,1.496)$, ${\mathbf{k}}_4=(1.496,1.496)$. Why? End slope.
5. Combine (both components identical): $y_{1,1}=1+\frac{0.4}{6}(1+2(1.2)+2(1.24)+1.496)=1+\frac{0.4}{6}(7.376)\approx \mathrm{1.491733.}$ Why this step? Weight the four positive slopes $1,2,2,1$ over $6$ and step forward; the doubled midpoints keep the exponential's accelerating curve accurate.

Verify: $e^{0.4}=1.491825$; error $\approx 9\times 10^{-5}$. Both slopes stayed positive as forecast.

---

Ex 7 — Cell G: limiting / stiff twist

Forecast: the true solution is $(e^{-20t},e^{-t})$; at $t=0.2$ that's $(e^{-4},e^{-0.2})=(0.0183,0.8187)$. But the fast component decays with rate $20$, and $h\cdot 20=4$ is large. RK4's stability limit for $y^{\prime }=\lambda y$ needs roughly $|h\lambda |\lesssim 2.78$. With $|h\lambda |=4$ we are outside it — expect the $y_1$ estimate to blow past reason (possibly wrong sign or magnitude), while $y_2$ (rate $1$, $|h\lambda |=0.2$) stays fine.

Compute $y_1$ (a scalar RK4 with $f=-20y_1$, $h=0.2$, so $h\lambda =-4$). For a linear $y^{\prime }=\lambda y$, one RK4 step multiplies by $R(z)=1+z+\frac{z^2}{2}+\frac{z^3}{6}+\frac{z^4}{24},z=h\lambda =-4.$

1. $z=-4$. Why? The step's stability argument.
2. $R(-4)=1-4+8-\frac{64}{6}+\frac{256}{24}=1-4+8-10.6667+10.6667=5.0$. Why? Plug into the polynomial.
3. So $y_{1,1}=1\cdot R(-4)=5.0$. Why? One step scales the state by $R(z)$ — this factor is the combined effect of the four-slope $1,2,2,1$ blend for a linear field.
4. For $y_2$: $z=-0.2$, $R(-0.2)=1-0.2+0.02-0.001333+0.0000667\approx 0.818733$. Why? Small $|z|$, safely stable.

Verify: true $y_1(0.2)=e^{-4}=0.0183$ but RK4 gave $5.0$ — a catastrophic overshoot (wrong by a factor $\sim 270$, and $|R(-4)|=5>1$ so it will grow each step). $y_2$ gave $0.8187$ vs true $e^{-0.2}=0.8187$ — perfect. Fix: shrink $h$ so $|h\lambda |<2.78$ for the stiffest rate, i.e. $h<2.78/20=0.139$, or use an implicit method. Explicit RK4 is limited by the fastest component.

---

Ex 8 — Cell H: non-autonomous (explicit $t$)

Forecast: now $\mathbf{f}$ mentions $t$ explicitly, so the $t_n+\frac{h}{2}$ and $t_n+h$ arguments actually matter — this is the case where forgetting to advance time breaks things. Start slope $y_2^{\prime }=0-0=0$, $y_1^{\prime }=1$: position rising, velocity flat. Because the exact answer is the straight line $y=t$, RK4 should return $0.2$ and $1$ with essentially no error.

Here $\mathbf{f}(t,\mathbf{y})=(y_2,t-y_1)$, state $(0,1)$ at $t_0=0$.

1. ${\mathbf{k}}_1=(1,0-0)=(1,0)$. Why? Slopes at $t=0$, state $(0,1)$.
2. $t_0+\frac{h}{2}=0.1$, state $=(0,1)+0.1(1,0)=(0.1,1)$; ${\mathbf{k}}_2=(1,0.1-0.1)=(1,0)$. Why? Midpoint — note we used $t=0.1$, not $0$.
3. State $=(0,1)+0.1(1,0)=(0.1,1)$; ${\mathbf{k}}_3=(1,0.1-0.1)=(1,0)$. Why? Refined midpoint at $t=0.1$.
4. $t_0+h=0.2$, state $=(0,1)+0.2(1,0)=(0.2,1)$; ${\mathbf{k}}_4=(1,0.2-0.2)=(1,0)$. Why? End, at $t=0.2$.
5. Combine: $y_{1,1}=0+\frac{0.2}{6}(1+2(1)+2(1)+1)=\frac{0.2}{6}(6)=\mathrm{0.2.}$ $y_{2,1}=1+\frac{0.2}{6}(0+0+0+0)=1.$ Why this step? The $1,2,2,1$ blend of the four slopes gives average velocity $(1,0)$, so we step by $h\cdot (1,0)=(0.2,0)$ — but only because we evaluated each ${\mathbf{k}}_i$ at the correct time.

Verify: the exact solution is $y(t)=t$ (check: $y^{\prime \prime }=0$, so $y^{\prime \prime }+y=0+t=t$ ✓; $y(0)=0$, $y^{\prime }(0)=1$ ✓), giving $y(0.2)=0.2$ and $y^{\prime }(0.2)=1$. RK4 returned exactly $y_{1,1}=0.2$ and $y_{2,1}=1$ — zero error. This is because we advanced $t$ correctly in each stage; had we frozen $t=0$, the $t-y_1$ term would be evaluated at the wrong time and the answer would be wrong.

---

Recall

Recall Quick self-test

Why did Ex 5 stay exactly put? ::: At the fixed point $\mathbf{f}=0$, so all four ${\mathbf{k}}_i$ vanish and the weighted average adds nothing. Why did Ex 7 return $5.0$ instead of near $0$? ::: The RK4 amplification $R(h\lambda )$ at $h\lambda =-4$ equals $5>1$, so the step amplifies instead of decaying — we were outside the stability region. In Ex 8, what would break if you kept $t=0$ in every stage? ::: The $t-y_1$ term would be evaluated at the wrong time, giving wrong $\mathbf{k}$'s and losing the exact answer.

---

Connections

• Parent: RK4 for systems
• RK4 for a single ODE — Ex 2 shows the decoupled case reduces to this.
• Reducing higher-order ODEs to first-order systems — Ex 3, Ex 8.
• Euler's method for systems — the cheaper cousin; try these examples with it to see the accuracy gap.
• Local vs Global Truncation Error — why bigger rates in Ex 2 gave bigger error.
• Stiff systems and stability — the failure in Ex 7.
• Simpson's Rule — the quadrature RK4 imitates.