Worked examples — Implementing ODE solvers from scratch — Euler, RK4

You already met Euler and RK4 in the parent note. This page does one thing: it drags those two update rules through every kind of situation a step-marching solver can face — growth, decay, oscillation, a curve that bends the "wrong" way, a step so big the method blows up, a real-world cooling problem, and an exam trap. If you can follow all of these, no ODE-solver question can surprise you.

Before we start, let me re-anchor the two tools in plain words so nothing below uses a symbol you haven't been handed.

The word "slope" is the hero of every example. Keep this picture in your head: a hiker in fog reading a compass at their feet (Euler) versus a hiker who peeks ahead three times before committing (RK4).

---

The scenario matrix

Every case class an ODE-stepper can be asked about, and the example that hits it:

Cell	What makes it distinct	Example
A. Growth $y^{\prime }=+ky$	slope grows, curve bends up → Euler undershoots	Ex 1
B. Decay $y^{\prime }=-ky$	slope shrinks, curve bends toward zero → sign of error flips	Ex 2
C. Non-autonomous $y^{\prime }=f(t)$	slope depends on $t$ — the $t$-shifts inside $k_2,k_3,k_4$ matter	Ex 3
D. Oscillation $y^{\prime }=\cos t$	curve turns both ways; slope passes through zero	Ex 4
E. Degenerate / zero slope $y^{\prime }=0$	flat line; both methods must return a constant	Ex 5
F. Big-step instability large $h$ on decay	both methods blow up past a critical $h$	Ex 6
G. Real-world word problem	Newton cooling, must set up $f$ and read units	Ex 7
H. Exam twist	reuse-$k_1$ trap + error-scaling reasoning	Ex 8

We work all eight below. Each one asks you to Forecast first — actually guess before reading the steps.

---

Example 1 — Cell A: exponential growth $y^{\prime }=y$, $y(0)=1$

Forecast: the curve bends upward (slope keeps increasing). Euler reads the slope only at the low-left start, so it uses too-small a slope → it should undershoot. Which method lands closer? Guess now.

1. Euler. $y_1=y_0+hf(t_0,y_0)=1+0.5\cdot (1)=1.5$. Why this step? Euler assumes the starting slope $f(0,1)=1$ holds for the whole half-step. Since the real slope rises during the step, $1.5$ falls short of $1.6487$.
2. RK4 $k_1$. $k_1=f(0,1)=1$. Why? Slope right where we stand — the same reading Euler used.
3. RK4 $k_2$. $k_2=f(0.25,1+0.25\cdot 1)=1.25$. Why? Peek at the midpoint $t=0.25$, nudging $y$ forward with $k_1$. The slope there is bigger — this already corrects the undershoot.
4. RK4 $k_3$. $k_3=f(0.25,1+0.25\cdot 1.25)=1.3125$. Why? Re-do the midpoint but nudge $y$ with the better slope $k_2$. Each stage refines the last.
5. RK4 $k_4$. $k_4=f(0.5,1+0.5\cdot 1.3125)=1.65625$. Why? The end-of-step slope, using $k_3$ to get there.
6. RK4 combine. $y_1=1+\frac{0.5}{6}(1+2(1.25)+2(1.3125)+1.65625)=1+\frac{0.5}{6}(7.71875)\approx 1.64323$. Why the $1,2,2,1$? Weighted mean slope; midpoints count double because they best represent the average over the step. Note the inner sum: $1+2.5+2.625+1.65625=7.78125$… wait — count carefully: $1+2(1.25)+2(1.3125)+1.65625=1+2.5+2.625+1.65625=7.78125$.

Let me state the arithmetic cleanly so nothing is fudged: the weighted sum is $k_1+2k_2+2k_3+k_4=1+2(1.25)+2(1.3125)+1.65625=7.78125,$ so $y_1=1+\frac{0.5}{6}(7.78125)=1+0.6484375=\mathrm{1.6484375.}$

Verify: true $e^{0.5}=1.64872$. Euler error $=0.1487$; RK4 error $=0.00028$ — RK4 is ~530× closer on a single step. Undershoot forecast confirmed for Euler. ✓

Read the figure: in Figure s01 the white curve is the true $e^t$. Follow the amber line — it is Euler's straight hop, and it ends below the white curve: that gap is the undershoot you forecast. The cyan square is RK4's endpoint, sitting practically on the white curve. The visual takeaway: one straight Euler segment cannot hug a bending curve, but RK4's four-slope average can.

---

Example 2 — Cell B: exponential decay $y^{\prime }=-y$, $y(0)=1$

Forecast: now the curve bends toward zero from above — it is convex, flattening as it drops. Euler uses the steep starting slope $-1$ for the whole step. Will it land above or below the true value?

1. Euler. $y_1=1+0.5\cdot (-1)=0.5$. Why this step? The starting slope $f(0,1)=-1$ is the steepest descent the curve ever has here; using it for the whole step drops us too far.
2. Compare. True $e^{-0.5}=0.6065$, so Euler at $0.5$ undershoots (lands below). Why note this? In Cell A (growth) Euler landed below the curve too, but there "below" meant "too small a rise". Here "below" means "dropped too fast". The lesson: Euler always errs on the side of the starting slope, and whether that's high or low depends on curvature — but for both growth and this decay it lands under the true curve because both are convex.

Verify: $|0.6065-0.5|=0.1065$ actual error. The leading-order Taylor estimate is $\frac{h^2}{2}|y^{\prime \prime }|=\frac{0.25}{2}\cdot 0.6065\approx 0.076$ — the right order of magnitude; higher-order terms in the expansion account for the gap up to the true $0.107$. ✓

---

Example 3 — Cell C: non-autonomous $y^{\prime }=2t$, $y(0)=0$

Forecast: here $f$ ignores $y$ entirely — the slope is set purely by time. This is exactly the case where forgetting the $t$-shifts inside $k_2,k_3,k_4$ silently ruins you. Guess: will RK4 be exact or approximate for a parabola?

Throughout this example the step size is $h=1$, so wherever you see $\frac{h}{2}$ it means $0.5$ and $\frac{h}{6}$ means $\frac{1}{6}$ — I keep $h$ in the formulas and substitute at the end so the general shape stays visible.

1. $k_1$. $k_1=f(0,0)=2\cdot 0=0$. Why? Slope at the very start, $t=0$, is zero.
2. $k_2$. $k_2=f(0+\frac{h}{2},0+\frac{h}{2}\cdot 0)=f(0.5,0)=2\cdot 0.5=1$. Why? We must evaluate at $t=0+\frac{h}{2}=0.5$, not $t=0$. The $y$-nudge is irrelevant because $f$ doesn't read $y$, but the $t$-shift is everything.
3. $k_3$. $k_3=f(0.5,0+\frac{h}{2}\cdot 1)=f(0.5,0.5)=2\cdot 0.5=1$. Why? Same $t=0.5$ → same slope $1$. Confirms the two midpoint slopes agree when $f$ is $y$-independent.
4. $k_4$. $k_4=f(0+h,\dots )=f(1,\dots )=2\cdot 1=2$. Why? End-of-step slope at $t=0+h=1$.
5. Combine. Keeping the factor $h$ explicit: $y_1=0+\frac{h}{6}(k_1+2k_2+2k_3+k_4)=\frac{1}{6}(0+2(1)+2(1)+2)=\frac{1}{6}(6)=1.$ Why? Weighted average of the four slopes, scaled by $\frac{h}{6}=\frac{1}{6}$ since $h=1$, lands exactly.

Verify: true $t^2$ at $t=1$ is $1$ → RK4 error $=0$. Why exact? RK4 matches the Taylor series through $h^4$; a quadratic has no terms beyond $h^2$, so RK4 is exact for it. ✓

---

Example 4 — Cell D: oscillation $y^{\prime }=\cos t$, $y(0)=0$

Forecast: the slope starts at $\cos 0=1$ (rising fast) but decreases to $0$ by the end of this quarter-period. Euler locks in the biggest slope of the interval. Overshoot or undershoot?

1. Euler. $y_1=0+\frac{\pi }{2}\cdot \cos 0=\frac{\pi }{2}\cdot 1\approx 1.5708$. Why this step? Euler trusts the start slope $1$ across the whole quarter-turn. But the real slope keeps shrinking to zero, so a constant slope of $1$ overshoots badly.
2. Interpret. True value $1$, Euler gives $1.5708$ → overshoot of $0.5708$. Why the flip from Ex 1? In growth the missed slope was larger later (undershoot). Here the missed slope is smaller later, so the same "trust-the-start" rule now overshoots. This is why curvature sign, not method, decides the direction of error.

Verify: $\frac{\pi }{2}-1=0.5708$. Overshoot confirmed. ✓

Read the figure: in Figure s02 the white curve is $\sin t$ rising to $1$ then flattening. The amber straight line is Euler's step — it climbs at the constant steep slope $1$ and shoots above the white curve, ending at $\pi /2\approx 1.57$ while the true endpoint (cyan square) sits at $1$. Contrast this with s01: same "trust the start slope" rule, but here the curve flattens so the error points the opposite way. That sign-flip is the whole lesson.

---

Example 5 — Cell E: degenerate zero slope $y^{\prime }=0$, $y(0)=7$

Forecast: slope is zero everywhere. A correct solver must return $7$ unchanged. If yours drifts, it's buggy.

1. Euler. $y_1=7+10\cdot 0=7$. Why? Zero slope × any step = no movement.
2. RK4. $k_1=k_2=k_3=k_4=0$, so $y_1=7+\frac{10}{6}(0)=7$. Why? Every peek reads zero; the weighted average of zeros is zero.

Verify: both return exactly $7$ for any $h$. ✓ This is the cheapest but most important unit test — run it first on any solver you write.

---

Example 6 — Cell F: big-step instability $y^{\prime }=-3y$, $y(0)=1$

First, a symbol we need. When an ODE has the form $y^{\prime }=\lambda y$, the number $\lambda$ (Greek "lambda") is just the growth/decay coefficient — the same role played by $\pm k$ in the matrix's rows "$y^{\prime }=+ky$" and "$y^{\prime }=-ky$". For our problem $y^{\prime }=-3y$ we simply read off $\lambda =-3$ (a decay, because it is negative). Everything below is written with $\lambda$ so the rule is reusable, but for us $\lambda =-3$.

Forecast: the true solution decays monotonically toward zero. But watch what a big step does. Guess: will Euler's numbers stay small and positive?

1. Step 1. $y_1=1+1\cdot (-3\cdot 1)=1-3=-2$. Why the disaster? Euler multiplies the step by the slope $-3$. With $h=1$ the correction $-3$ overshoots past zero into negative territory — physically impossible for a decay.
2. Step 2. $y_2=-2+1\cdot (-3\cdot -2)=-2+6=4$. Why? The overshoot flips the sign, the slope flips too, and the error grows each step.
3. Step 3. $y_3=4+1\cdot (-3\cdot 4)=4-12=-8$. Why? Each step the magnitude multiplies by $|1+h\lambda |=|1+1\cdot (-3)|=|-2|=2$. It doubles forever — unstable.

Euler's stability rule: forward Euler on $y^{\prime }=\lambda y$ stays bounded only if $|1+h\lambda |\le 1$. Here $\lambda =-3$, so we need $|1-3h|\le 1\Rightarrow 0\le h\le \frac{2}{3}$. Our $h=1$ violates it.

Does RK4 escape? No — higher order does not mean unconditionally stable. For $y^{\prime }=\lambda y$ each RK4 step multiplies $y$ by the polynomial $R(z)=1+z+\frac{z^2}{2}+\frac{z^3}{6}+\frac{z^4}{24},z=h\lambda ,$ (the Taylor series of $e^z$ truncated at $z^4$ — exactly the four stages). RK4 is bounded only if $|R(z)|\le 1$. With $h=1,\lambda =-3$ we get $z=-3$ and $R(-3)=1-3+4.5-4.5+3.375=1.375,$ whose magnitude $>1$, so RK4 with $h=1$ also blows up on this problem. RK4's stability limit here is roughly $h\lesssim 0.93$ (where $|R(h\lambda )|=1$) — larger than Euler's $\frac{2}{3}$, but still finite. The lesson: for stiff decay you must shrink $h$ for both methods; see Numerical Stability and Stiff ODEs.

Verify: with the safe $h=0.5$: Euler $y_1=1+0.5(-3)=-0.5$ — still negative but $|1-1.5|=0.5\le 1$, so magnitudes now shrink, not grow. Sequence $-0.5,0.25,-0.125,\cdots \to 0$. Bounded. For RK4 at $h=0.5$, $z=-1.5$ gives $R(-1.5)=1-1.5+1.125-0.5625+0.2109=0.2734$, $|R|<1$, so RK4 is bounded too. ✓

Read the figure: in Figure s03 the white curve is the true decay $e^{-3t}$ hugging zero. The amber zig-zag is Euler with $h=1$: it flings above and below the white dashed zero-line, each swing bigger than the last — that is the doubling you computed. The cyan markers are Euler with the safe $h=0.5$: same alternating signs but each swing smaller, spiralling into zero. The picture makes "unstable vs bounded" a matter of whether the zig-zag grows or shrinks.

---

Example 7 — Cell G: real-world word problem (Newton's cooling)

Forecast: the cup cools fast at first (big gap → steep slope), then slowly. Over a whole 10-minute step, Euler will use only the initial steep slope. Overshoot or undershoot the true temperature?

Set up: $f(t,T)=-0.1(T-20)$. Units: $[{}^{\circ }\text{C}/\text{min}]$. Start $T_0=90$.

1. Euler. slope $=-0.1(90-20)=-7{}^{\circ }\text{C/min}$. $T_1=90+10\cdot (-7)=90-70=20$. Why suspicious? Euler predicts the coffee reaches room temperature exactly, using the fast initial cooling for all 10 minutes. Real cooling slows down, so it can't cool that much.
2. RK4 $k_1$. $k_1=-0.1(90-20)=-7$.
3. RK4 $k_2$. $T$-nudge: $90+5\cdot (-7)=55$; $k_2=-0.1(55-20)=-3.5$. Why? midpoint slope is gentler because the cup already cooled.
4. RK4 $k_3$. $90+5\cdot (-3.5)=72.5$; $k_3=-0.1(72.5-20)=-5.25$. Why? refine the midpoint with $k_2$.
5. RK4 $k_4$. $90+10\cdot (-5.25)=37.5$; $k_4=-0.1(37.5-20)=-1.75$. Why? end slope, gentlest of all.
6. Combine. The weighted sum of slopes is $k_1+2k_2+2k_3+k_4=-7+2(-3.5)+2(-5.25)+(-1.75)=-7-7-10.5-1.75=-26.25$. So $T_1=90+\frac{10}{6}(-26.25)=90-43.75=46.25{}^{\circ }\text{C}$.

Verify: exact solution $T(t)=20+70e^{-0.1t}$, so $T(10)=20+70e^{-1}\approx 20+25.7517=45.7517{}^{\circ }\text{C}$. Euler's $20^{\circ }$ is a gross overshoot in cooling (off by $25.75^{\circ }$); RK4's $46.25^{\circ }$ is off by only $0.498^{\circ }$ despite one giant 10-minute step. ✓ Units all in ${}^{\circ }\text{C}$. Forecast (Euler cools too much) confirmed.

---

Example 8 — Cell H: exam twist (the reuse-$k_1$ trap + scaling)

Forecast (a): reusing $k_1$ makes $k_2$ and $k_3$ identical — the method loses an independent slope. Will $y_1$ be too big or too small versus the correct $2.7083$?

Part (a):

1. Correct stages (from parent): $k_1=1,k_2=1.5,k_3=1.75,k_4=2.75$, giving correct $y_1=1+\frac{1}{6}(1+3+3.5+2.75)=1+\frac{10.25}{6}\approx 2.70833$.
2. Buggy $k_3$. $k_3=f(0.5,1+0.5\cdot k_1)=f(0.5,1.5)=1.5$ (same as $k_2$). Why the damage? $k_3$ no longer refines $k_2$; the two midpoint peeks now carry identical information, so the averaging is weaker.
3. Buggy $k_4$. $k_4=f(1,1+1\cdot k_3)=f(1,2.5)=2.5$. Why changed? $k_4$ builds on the wrong $k_3$.
4. Combine. $y_1=1+\frac{1}{6}(1+2(1.5)+2(1.5)+2.5)=1+\frac{1}{6}(9.5)\approx 2.58333$. Why worse? True $e=2.71828$; correct RK4 gave $2.70833$ (error $0.00995$), buggy gives $2.58333$ (error $0.135$) — the bug is ~$13.5\times$ worse, collapsing toward Euler-like quality.

Part (b):

1. Scaling law. RK4 error $\propto h^4$. Going from $h=0.1$ to $h=0.025$ is a factor $\frac{0.1}{0.025}=4$ reduction in step. Why $h^4$? Global error order is $O(h^4)$ (parent), so error shrinks by $4^4=256$.
2. Predict. $2.1\times 10^{-6}/256\approx 8.2\times 10^{-9}$.

Verify: (a) buggy inner sum $\frac{9.5}{6}=1.58333$, so $y_1=2.58333$; its gap from the correct RK4 value $2.70833$ is $0.125$. ✓ (b) $2.1\times 10^{-6}/256=8.20\times 10^{-9}$. ✓

Recall Which cell taught which lesson?

Growth vs decay :::: curvature sign, not the method, decides over/undershoot direction (Ex 1, 2, 4). Non-autonomous :::: never skip the $t$-shift inside $k_2,k_3,k_4$ (Ex 3). Zero slope :::: both methods must return a constant for any $h$ — your first unit test (Ex 5). Big step :::: Euler is bounded only if $|1+h\lambda |\le 1$; RK4 is bounded only if $|R(h\lambda )|\le 1$ — both have a finite limit (Ex 6). Reuse-$k_1$ :::: each stage must feed the next or you drop back toward first order (Ex 8a).

---

Connections

• Taylor Series Expansion — why RK4 is exact on the quadratic in Ex 3, and the source of the $\frac{h^2}{2}{y}^{\prime \prime }$ error in Ex 2.
• Simpson's Rule — the $1,2,2,1$ averaging philosophy behind Ex 1's RK4 step.
• Finite Difference Approximation of Derivatives — the forward difference that is Euler.
• Numerical Stability and Stiff ODEs — the blow-up in Ex 6 explained in full.
• Adaptive Step Size (RK45 / Dormand–Prince) — how solvers pick $h$ so Ex 6 never happens.
• scipy.integrate.solve_ivp — the production tool that wraps all of this.