Ordinary Differential Equations

From a law of change to a predicted future — the mathematics every engineer uses daily, whether or not they name it.

You almost never get a formula for where a system will be. You get a law for how it responds right now: a hot drink cools in proportion to how much warmer it is than the room; a spring pushes in proportion to how far it is stretched; a capacitor's voltage drifts at a rate set by the current through the resistor. An ordinary differential equation (ODE) is the device that turns such a "law of change" into a full trajectory — a prediction of the future from the present state. This page builds that idea in four moves: what the law says, what it looks like geometrically, how to solve the canonical cases exactly, and what to do when no formula exists (which is most of the time).

1The problem an ODE solves

Say you pour coffee at $T_0 = 95\,^\circ\mathrm{C}$ into a room held at $T_{\mathrm{env}} = 20\,^\circ\mathrm{C}$. Nobody hands you $T(t)$, the temperature as a function of time $t$. What you can state — from Newton's law of cooling — is how the temperature responds at this instant: the cooling rate is proportional to the current gap between the coffee and the room. In symbols, with $T(t)$ the temperature in $^\circ\mathrm{C}$, $t$ the time in minutes, and $k>0$ a constant measuring how vigorously heat escapes (units $\mathrm{min}^{-1}$):

$$\frac{dT}{dt} \;=\; -k\,\bigl(T - T_{\mathrm{env}}\bigr)$$

Why this form is universal: physical laws are local in time — causality means the future is driven by the present state, not by the future itself. So laws naturally relate a system's state and its instantaneous rate of change. Written as an equation for the derivative, that relation is an ODE: ordinary because the unknown $T$ depends on a single variable $t$ (no partial derivatives), differential because derivatives of the unknown appear, and an equation because it must hold at every instant.

Solving the ODE means inverting the relationship: from a rule about the derivative, recover the function $T(t)$ itself. You already do this instinctively — "distance = integral of velocity" is the same move. An ODE just makes it systematic when the velocity itself depends on where you are.

Takeaway. An ODE is a law of change — a rule for the derivative. Solving it turns that rule into a trajectory, and the initial state turns the trajectory into a prediction.

2Vocabulary: classify before you solve

Why classify first: the type of an ODE dictates both the method of attack and the behavior you should expect. Four distinctions carry most of the weight.

Unknown and variable. The unknown is a function — call it $y(x)$ — and $x$ is the independent variable (often time $t$). Notation: $y' = dy/dx$, $y'' = d^2y/dx^2$, and $\dot{y} = dy/dt$ when the variable is time.
Order. The order is the highest derivative that appears. $y' = f(x,y)$ is first order; $y'' + 2y' + y = 0$ is second order. Why it matters: an $n$-th order ODE needs $n$ initial conditions to pin down a unique solution — a spring needs both position and velocity to know what happens next.
Linear vs. nonlinear. The equation is linear if $y$ and its derivatives appear only to the first power and never inside other functions — the general first-order form is $y' + p(x)\,y = q(x)$. Anything with $y^2$, $\sin y$, or $y\,y'$ is nonlinear. Why it matters: linear equations admit superposition (combine solutions) and a complete toolbox; nonlinear ones have few formulas but rich behavior — equilibria, oscillation, chaos — and are handled geometrically and numerically.
Initial value problem (IVP). The ODE plus starting data, $y(x_0) = y_0$. Why needed: the law alone specifies a whole family of trajectories; the initial condition selects one, just as a navigation app needs your current position before it can route you.

$$y' = f(x,\,y)$$ one rule for the derivative — infinitely many solution curves

Linear

$y$ appears only to the first power; superposition works.

$y' = k\,y$ — growth & decay §4

$y' + p(x)\,y = q(x)$ — integrating factor §4.3

$y'' + 2\zeta\omega_n y' + \omega_n^2\,y = 0$ — vibration §6

driven by forcing $F(t)$ — resonance §7

Nonlinear

Products and functions of $y$; few formulas, rich behavior.

$y' = y(1-y)$ — separable, logistic §5

$y' = y^2$ — finite-time blow-up §3

everything else — read the field, integrate numerically §8

A map of this page: every ODE you meet is first sorted into one of these bins, and each bin links to the section that handles it.

Takeaway. Order tells you how many initial conditions you owe; linearity tells you whether algebra (superposition) or geometry (field reading) is the sharper tool.

3The geometric picture: an ODE is a field of instructions

Why look at geometry before hunting formulas: a formula answers "what is $y$ at $x$?", but the ODE itself answers a more local question — "if I am at state $(x, y)$, what must my slope be?" Take the concrete equation

$$y' = x - y.$$

At every point $(x,y)$ of the plane, the right-hand side $x - y$ assigns one number: the slope any passing solution must have there. Draw a short segment with that slope at each point of a grid and the plane becomes a direction field — a set of local instructions, like wind arrows on a weather map. A solution is any curve that is tangent to the instruction at every point it crosses.

Direction field of y' = x - y with three solution curves — **Figure 1 — The direction field of $y' = x - y$.** Gray segments are the local instructions (slope $= x - y$ at each point). Starting from the three black dots $y(0) = 2.5,\ 0,\ -0.8$ and always following the segments produces the three colored solution curves. All of them are funneled onto the dashed line $y = x - 1$, where the instruction is $y' = 0$ — a *moving equilibrium* that every trajectory approaches.

Three consequences follow directly from this picture, and each is a load-bearing fact:

An IVP is "start here." Pick the starting point; the field determines the rest of the curve. Why you may trust this: the existence and uniqueness theorem (Picard–Lindelöf) makes it precise — if $f(x,y)$ is continuous and does not change too abruptly in the $y$ direction (its $y$-slope is bounded, the Lipschitz condition), then exactly one solution passes through each starting point, at least locally. Why engineers care: this is the license for simulation — one measured state determines the trajectory, and small measurement errors produce small trajectory errors (the problem is well-posed).
Uniqueness is local, existence may be finite. For $y' = y^2$ the solution $y = 1/(c - x)$ escapes to infinity in finite time. Nonlinear fields can eject solutions — the red curve in Figure 6 does exactly this. "The solution exists" is a theorem, not a given.
Equilibria are visible. Wherever the instruction is horizontal ($y' = 0$), a solution may rest forever. Finding and classifying these rest points — before solving anything — is often the whole engineering question (will it settle? diverge? oscillate?).

Why this picture also sets up numerics: "follow the segments" is an algorithm waiting to happen — Section 8 automates exactly this threading.

Takeaway. An ODE turns the plane into a slope field; a solution is a curve threading it, and the initial condition says where to start threading. Smooth fields give one well-defined thread per starting point.

4First-order linear: the exponential machinery

Why this section earns its keep: "rate proportional to amount" is arguably the most common law in engineering — cooling, radioactive decay, RC discharge, population growth, compound interest. Master this one case and a large fraction of practice is covered.

4.1 The pure case $y' = k\,y$

Why expect an exponential: the exponential is the function that reproduces itself under differentiation — $\frac{d}{dx}e^{kx} = k\,e^{kx}$. A law saying "the derivative is proportional to the function" is a job advertisement for exactly that function. Rather than guess, derive it, with a move worth internalizing:

Step 1 — divide out the growth. Multiply the equation by $e^{-kx}$ and use the product rule:

$$\frac{d}{dx}\Bigl(y\,e^{-kx}\Bigr) \;=\; e^{-kx}\,\underbrace{\bigl(y' - k\,y\bigr)}_{=\;0\ \text{by the ODE}} \;=\; 0.$$

Why this works: the factor $e^{-kx}$ was chosen precisely so that the leftover bracket $y' - ky$ is the ODE's left-hand side — the multiplier converts the whole equation into "some derivative is zero."

Step 2 — a function with zero derivative is constant. That is the fundamental theorem of calculus read backwards, so $y\,e^{-kx} = C$, hence

$$y(x) = C\,e^{kx}.$$

Step 3 — spend the initial condition. $y(0) = y_0$ forces $C = y_0$, so the IVP's answer is $y = y_0\,e^{kx}$. The family $C\,e^{kx}$ is the one-parameter story the ODE tells; the IC picks the member.

Read the constant $k$ as a relative growth rate: $y'/y = k$ says each unit of $y$ grows at the same fractional pace. Its reciprocal $\tau = 1/k$ is the time constant — the time for the amount to change by a factor of $e \approx 2.718$. Negative $k$ is decay with the same clock.

4.2 Worked example: cooling coffee, start to numbers

Back to Section 1's coffee: $T' = -k\,(T - 20)$, $T(0) = 95$. Why shift variables first: the room temperature is a constant offset; measuring the gap makes the law pure. Let $G = T - 20$ (the gap, in $^\circ\mathrm{C}$). Then $G' = T'$ (constants don't differentiate), so $G' = -k\,G$ — Section 4.1 verbatim:

$$G(t) = 75\,e^{-kt} \qquad\Longrightarrow\qquad T(t) = 20 + 75\,e^{-kt}.$$

Calibrate $k$ from one measurement. Suppose the coffee reads $65.5\,^\circ\mathrm{C}$ after $t = 5$ min. Then $G(5) = 45.5 = 75\,e^{-5k}$, and solving for $k$ uses the logarithm — the inverse of the exponential:

$$k \;=\; \frac{1}{5}\,\ln\!\frac{75}{45.5} \;\approx\; 0.100\ \mathrm{min}^{-1}.$$

Predict. With $k = 0.1$: at $t = 10$ min, $T = 20 + 75 e^{-1} \approx 47.6\,^\circ\mathrm{C}$; the coffee cools to $30\,^\circ\mathrm{C}$ at $t = \ln(75/10)/0.1 \approx 20.1$ min. Note the time constant $\tau = 1/k = 10$ min doing its job: every 10 minutes the gap shrinks by a factor of $e$ — $75 \to 27.6$ at $t=10$, exactly what the formula gives.

Cooling curves on a linear scale — **Figure 2 — Exponential decay, two views.** (a) All curves start at $95\,^\circ\mathrm{C}$ and sag toward the ambient line; larger $k$ sags faster. (b) The same data with the gap on a log scale: every curve becomes a straight line whose slope is $-k$. Engineering habit worth keeping: exponential behavior betrays itself as a straight line on semi-log paper — use it to spot the law and read off $k$ from two data points.

The cooling gap on a logarithmic scale — straight lines — **Figure 2 — Exponential decay, two views.** (a) All curves start at $95\,^\circ\mathrm{C}$ and sag toward the ambient line; larger $k$ sags faster. (b) The same data with the gap on a log scale: every curve becomes a straight line whose slope is $-k$. Engineering habit worth keeping: exponential behavior betrays itself as a straight line on semi-log paper — use it to spot the law and read off $k$ from two data points.

The same law with different letters runs half of electronics: a discharging RC circuit obeys $V_c' = -V_c/(RC)$ for the capacitor voltage $V_c$, i.e. the pure case with $k = 1/RC$ — the familiar $\tau = RC$ time constant falls out as $1/k$.

4.3 The general linear ODE and the integrating factor

Why generalize: real systems have inputs and offsets — a heater also runs while the coffee cools, forcing terms appear. The general first-order linear ODE is

$$y' + p(x)\,y = q(x),$$

with $p, q$ known functions of $x$ alone. Why a special trick is needed: if $q \neq 0$ the pure-case move fails — $y'$ and $y$ no longer share a common factor. So we manufacture one: find a multiplier $\mu(x)$, called the integrating factor, such that the left side becomes a single product-rule derivative:

$$\mu\,y' + \mu\,p\,y \;\overset{!}{=}\; (\mu\,y)' \;=\; \mu\,y' + \mu'\,y \qquad\Longleftrightarrow\qquad \frac{\mu'}{\mu} = p(x).$$

That last equation is itself the pure case ($\mu' = p\,\mu$) — the trick eats its own dog food — so $\mu = e^{\int p\,dx}$ (any antiderivative; a constant factor would cancel). The ODE collapses to one integration:

$$(\mu\,y)' = \mu\,q \qquad\Longrightarrow\qquad y = \frac{1}{\mu}\left(\int \mu\,q\;dx + C\right).$$

Worked example — closing the loop with Figure 1. The direction-field equation $y' = x - y$ rearranges to $y' + y = x$, so $p = 1$, $q = x$, $\mu = e^{x}$. Then $(y\,e^x)' = x\,e^x$; integrating by parts ($\int x e^x dx = x e^x - \int e^x dx$, because $\frac{d}{dx}(x-1)e^x = x e^x$) gives $y\,e^x = (x-1)e^x + C$, so

$$y = x - 1 + C\,e^{-x}.$$

The blue curve in Figure 1 started at $y(0) = 2.5$, i.e. $C = 3.5$: $y = x - 1 + 3.5\,e^{-x}$. Formula and field agree — and the formula explains the attractor seen in the figure: as $x$ grows, the term $C\,e^{-x}$ dies, and every solution collapses onto the line $y = x - 1$.

Takeaway. Rate proportional to amount gives an exponential; the integrating factor manufactures a product-rule derivative so that any first-order linear ODE becomes a single integration.

5Nonlinear but separable: the logistic law and phase-line thinking

Why leave linear land: exponential growth cannot last — any real population, market, or reaction runs into a resource limit. The logistic equation adds the missing feedback: growth slows as the room runs out,

$$y' \;=\; r\,y\left(1 - \frac{y}{K}\right),$$

with $r > 0$ the intrinsic growth rate and $K > 0$ the carrying capacity. It is nonlinear (the factor $y(1-y/K)$ is quadratic in $y$) — but it is separable: the right side factors into a piece of $x$ (here, the constant $r$) times a piece of $y$.

Why separation is legal: for $y' = g(x)\,h(y)$, define $H(y) = \int dy/h(y)$. By the chain rule, $\frac{d}{dx}H\bigl(y(x)\bigr) = H'(y)\,y' = y'/h(y) = g(x)$ — so integrating both sides in $x$ is exactly the informal "collect the $y$'s with $dy$, the $x$'s with $dx$." The dance of differentials is bookkeeping for the chain rule.

Worked example — solve the logistic equation. Separating and using partial fractions (checkable by putting the pieces over the common denominator $y(1-y/K)$):

$$\int \frac{dy}{y\,(1-y/K)} \;=\; \int r\,dx \qquad\Longrightarrow\qquad \ln\left|\frac{y}{1-y/K}\right| = r\,x + c \qquad\Longrightarrow\qquad y(x) = \frac{K}{1 + A\,e^{-rx}},$$

where the constant $A = (K - y_0)/y_0$ encodes the initial condition. Why the S-curve shape matters: early on ($y \ll K$) the $1 - y/K$ factor is $\approx 1$ and growth is exponential; late ($y \to K$) the gap to capacity drives an exponential approach — the two pure cases of §4 glued together. Concrete timing: seeded at $y_0 = 0.1\,K$, the population reaches $90\%$ of capacity when $1 + 9e^{-rx} = 10$, i.e. after $x = \ln 81 / r \approx 4.39/r$.

Logistic direction field, S-curves from several initial conditions, equilibria at 0 and K — **Figure 6 — The logistic field $y' = y(1-y)$ (here $r = K = 1$).** The dashed lines mark the equilibria $y = 0$ and $y = 1$ where the instructions go horizontal. Solutions seeded between 0 and 1 (blue, orange) ride the S-curve up to $K$; seeded above $K$ (green) they sink back to it. Seeded below zero (red), the growth term $y(1-y)$ is negative and the solution dives to $-\infty$ in finite time — the denominator $1 + A e^{-rx}$ hits zero when $A < 0$. Same equation, three fates, decided entirely by which side of the unstable equilibrium you start on.

Phase-line reading — the skill that survives unsolvable equations. Why practice this: the exact formula above is pleasant, but the qualitative conclusions need no integral at all, and they generalize to equations you cannot solve. For an autonomous law $y' = f(y)$: mark the roots of $f$ (equilibria); check the sign of $f$ between roots (the direction of motion); equilibria where the sign flips from $+$ to $-$ are stable (arrows converge, like $y = K$); from $-$ to $+$ they are unstable (arrows diverge, like $y = 0$). The logistic field, read this way, explains all three fates in Figure 6 in ten seconds.

Beyond biology, the same equation models technology adoption, autocatalytic reactions, and epidemic early growth — and the sign-reading habit is the doorway to phase portraits and stability analysis in control theory.

Takeaway. Separable equations yield to two integrals; the deeper prize is reading equilibria and their stability straight off the sign of $f(y)$ — no formula required.

6Second order: vibration, damping, and the characteristic equation

Why second order deserves its own section: Newton's $F = ma$ is a statement about acceleration — the second derivative — so mechanics speaks second-order ODEs. The atom of vibration is the mass–spring–damper:

$$m\,y'' + c\,y' + k\,y = 0,$$

where $y(t)$ is displacement from rest, $m$ the mass, $k$ the spring stiffness (restoring force $-k\,y$), and $c$ the damping coefficient (dissipative force $-c\,y'$). Dividing by $m$ and renaming constants gives the engineer's normalized form

$$y'' + 2\zeta\,\omega_n\,y' + \omega_n^2\,y = 0,\qquad \omega_n = \sqrt{k/m},\quad \zeta = \frac{c}{2\sqrt{m\,k}},$$

with $\omega_n$ the natural frequency (rad/s) and $\zeta$ the dimensionless damping ratio. $\zeta = 0$ means undamped, $\zeta = 1$ is called critical damping. (The figures below take $\omega_n = 1$, $y(0) = 1$, $y'(0) = 0$: displace the mass and let go.)

Why try $y = e^{rt}$: with constant coefficients, differentiation acts on $e^{rt}$ as multiplication by $r$: $\frac{d}{dt}e^{rt} = r\,e^{rt}$. Each derivative just tags on another factor $r$, so the ODE becomes algebra:

$$\bigl(r^2 + 2\zeta\omega_n\,r + \omega_n^2\bigr)\,e^{rt} = 0 \qquad\Longrightarrow\qquad r^2 + 2\zeta\omega_n\,r + \omega_n^2 = 0.$$

This is the characteristic equation. Since $e^{rt}$ never vanishes, it forces the polynomial to vanish — calculus has been converted into root-finding. Why two roots suffice: a second-order IVP needs two initial conditions, its solution space is two-dimensional, and each distinct root contributes one independent exponential; the general solution is their sum (a two-dimensional vector space needs two basis vectors). Why complex roots still give real motion: Euler's formula $e^{i\theta} = \cos\theta + i\sin\theta$ converts complex exponentials into sines and cosines; taking real parts (legal by linearity) yields real solutions.

The quadratic formula gives $r = \omega_n\bigl(-\zeta \pm \sqrt{\zeta^2 - 1}\bigr)$ — the discriminant sorts all behavior into three regimes:

Regime	Roots $r$	Motion
$\zeta < 1$ underdamped	$-\zeta\omega_n \pm i\,\omega_d$, $\omega_d = \omega_n\sqrt{1-\zeta^2}$	decaying oscillation $e^{-\zeta\omega_n t}(A\cos\omega_d t + B\sin\omega_d t)$
$\zeta = 1$ critically damped	double root $-\omega_n$	$(A + B\,t)\,e^{-\omega_n t}$ — fastest return without overshoot
$\zeta > 1$ overdamped	two negative reals	slow exponential + fast exponential

Three damping regimes for y'' + 2 zeta y' + y = 0 — **Figure 3 — One mass–spring–damper, three personalities** ($y'' + 2\zeta y' + y = 0$, released from $y=1$ at rest). Dashed gray: $\zeta = 0$ oscillates forever — damping is the only thing that kills ringing. Blue ($\zeta = 0.15$): rings at $\omega_d \approx 0.99$ inside an $e^{-0.15t}$ envelope. Orange ($\zeta = 1$): the critically damped slide home — the fastest approach *with no overshoot*, which is why door closers and instrument needles are tuned near this value. Green ($\zeta = 2.5$): overdamped — its fast mode $e^{-4.79t}$ dies almost instantly, leaving the sluggish slow mode $e^{-0.21t}$ to crawl back; more damping, slower recovery.

Worked mini-example — do it once by hand. Solve $y'' - 3y' + 2y = 0$ with $y(0) = 1,\ y'(0) = 0$. Characteristic equation: $r^2 - 3r + 2 = (r-1)(r-2) = 0$, roots $r = 1, 2$, so $y = A e^{t} + B e^{2t}$. The initial conditions give $A + B = 1$ and $A + 2B = 0$, hence $B = -1,\ A = 2$:

$$y = 2e^{t} - e^{2t}.$$

Why always finish with a substitution check: differentiating, $y' = 2e^t - 2e^{2t}$, $y'' = 2e^t - 4e^{2t}$, and $y'' - 3y' + 2y = (2-6+2)e^{t} + (-4+6-2)e^{2t} = 0$ — plus $y(0) = 2 - 1 = 1$ and $y'(0) = 0$. Ten seconds of algebra buys certainty that the algebra above never slipped.

Takeaway. Constant coefficients turn calculus into algebra: exponentials are eigenfunctions of $d/dt$, and the roots of the characteristic polynomial are the decay rates and oscillation frequencies of the system.

7Forcing and resonance: the driven oscillator

Why add a drive: real systems are pushed — ground shakes a building, a signal drives a circuit, wind loads a bridge. Attach an external force to the mass–spring–damper:

$$y'' + 2\zeta\,\omega_n\,y' + \omega_n^2\,y = \frac{F_0}{m}\cos(\omega\,t).$$

Superposition splits the problem. Why legal: linearity — if $y_1$ solves the equation with zero right side and $y_p$ solves it with the drive, then substituting $y_h + y_p$ shows the cross-terms cancel and each piece does its own job. So the general solution is

$$y(t) \;=\; \underbrace{y_h(t)}_{\text{transient: dies out for } \zeta>0} \;+\; \underbrace{y_p(t)}_{\text{steady state: tracks the drive}}.$$

Long after startup, what remains is the steady response. Why guess $y_p = M\cos(\omega t - \varphi)$: a damped linear system driven at a single frequency responds at that same frequency, merely shifted in phase and scaled — plausible physically (the system cannot invent new frequencies) and confirmed by substitution; using $\cos = \operatorname{Re}(e^{i\omega t})$ turns the substitution into complex algebra. Carrying it out gives the amplitude ratio — output amplitude over the static deflection $F_0/k$ you would get by pushing slowly:

$$|H(i\omega)| \;=\; \frac{1}{\sqrt{\bigl(1 - (\omega/\omega_n)^2\bigr)^2 + \bigl(2\zeta\,\omega/\omega_n\bigr)^2}}.$$

Resonance curves: amplitude ratio vs frequency ratio for several damping ratios — **Figure 4 — Resonance curves of the driven oscillator** (log scale, as in a Bode magnitude plot). At low frequency all curves sit at 1: the mass follows the force quasi-statically. Near $\omega/\omega_n = 1$ the light-damping curves peak — for $\zeta = 0.05$ the amplification is $\approx 10\times$ at $\omega_r \approx 0.997\,\omega_n$. Above the peak, all curves fall in parallel — the $-40\,\mathrm{dB/decade}$ roll-off (mass-controlled region) you may recognize from filter design. Higher damping lowers and broadens the peak; $\zeta = 1$ does not peak at all.

Why the peak exists: near the natural frequency the drive pushes in phase with the mass's velocity each cycle, so it keeps adding energy; damping removes energy at a rate that grows with amplitude, and the balance sets a finite peak. Maximizing $|H|$ gives the exact resonance values (valid for $\zeta < 1/\sqrt{2}$):

$$\omega_r = \omega_n\sqrt{1 - 2\zeta^2}, \qquad |H|_{\max} = \frac{1}{2\zeta\sqrt{1 - \zeta^2}} \;\approx\; \frac{1}{2\zeta} \;=\; Q.$$

The approximation $1/2\zeta$ for light damping is the famous quality factor $Q$: damping ratio 1% means a $50\times$ amplification lurking at the natural frequency. This is why rotating machinery is balanced, why structures are modally tested before commissioning, and why the same mathematics gives bandwidth and overshoot in control design. Where this fits in the larger toolbox: $|H(i\omega)|$ is the magnitude of the transfer function of the second-order system; the Laplace transform (differentiation $\leftrightarrow$ multiplication by $s$) is the machinery that produces $H(s)$ from the ODE in one line — the algebraic twin of everything in §6.

Takeaway. Linear plus driven splits into a dying transient plus a steady tracking motion; the steady amplitude is the transfer function, and its peak — resonance — marks where energy injection and dissipation balance.

8When no formula exists: numerical integration

Why numerics is not a footnote: most ODEs — the three-body problem, almost every nonlinear circuit, real chemistry kinetics — have no closed-form solution at all. The field picture of §3 tells you what to do anyway: follow the instructions in small steps. That is precisely what a numerical solver is.

Euler's method — the simplest possible scheme. Stand at $(x_n, y_n)$, read the slope $f(x_n, y_n)$, walk a small step $h$ along it, and re-read:

$$y_{n+1} \;=\; y_n + h\,f(x_n,\,y_n).$$

Why justified: it is the definition of the derivative frozen at one point — the first-order Taylor approximation. Each step follows a straight segment while the true solution bends.

Worked example — $y' = y$, $y(0) = 1$, step $h = 0.5$ (Figure 5, orange): $y_1 = 1 + 0.5\cdot 1 = 1.5$; $y_2 = 1.5 + 0.5\cdot 1.5 = 2.25$; and so on to $y_4 = 1.5^4 = 5.0625$ at $x = 2$, versus the exact $e^2 \approx 7.389$. Every step uses the slope from the left edge of the interval, and since $e^x$ is convex, every straight chord undercuts the curve — Euler always lags here.

Euler's method steps versus the exact solution of y' = y — **Figure 5 — Euler's method is threading the direction field by hand.** Each colored segment walks the slope instruction at its left endpoint; dots mark the re-reads. Halving the step from $h = 0.5$ (orange) to $h = 0.25$ (green) visibly hugs the exact curve $y = e^x$ (dark) — and cuts the final error roughly in half, from $2.33$ to $1.43$.

How the error scales — count it. One step's error is $O(h^2)$ (the curvature term that Taylor's first order drops); a run to fixed $x$ takes $1/h$ steps, so the accumulated error is $O(h)$. The worked table confirms the half-when-$h$-halves law:

step $h$	Euler's $y(2)$ for $y' = y$	error vs. $e^2 \approx 7.3891$
0.5	5.0625	2.33
0.25	5.9605	1.43
0.125	6.5833	0.81
0.0625	6.9587	0.43

Why you will never actually use Euler: each refinement buys accuracy more cheaply than brute force:

Runge–Kutta 4 (RK4) samples the field four times per step (start, two midpoints, end) and combines them so the $h^2$ and $h^3$ error terms cancel — global error $O(h^4)$: roughly a million-fold accuracy gain for four times the work. This is the engine inside ode45 and scipy.integrate.solve_ivp.
Adaptive steps estimate the local error each step and shrink $h$ only where the field bends hard — no human babysitting the resolution.
Stiff problems mix fast and slow time constants (chemical kinetics, circuits with ns and ms modes). Explicit methods are hostage to the fastest mode — stability requires $h \lesssim 2/|\lambda_{\max}|$ — even when nothing interesting happens at that timescale. Implicit methods (e.g. backward Euler, $y_{n+1} = y_n + h\,f(x_{n+1}, y_{n+1})$, which solves a small equation per step) are unconditionally stable for decay modes, since their amplification factor $1/(1 + h\lambda)$ stays below 1 for any $h > 0$. Why this matters: picking explicit vs. implicit is often the single biggest performance decision in simulation.

Takeaway. A numerical solver is the slope-field picture turned into an algorithm; the step size is a resolution dial, and the error order tells you what each turn of the dial buys.

9Where this sits in the larger map

Why zoom out: each idea above is a door into a larger room; knowing the address makes the next steps findable.

Systems and state space. Real devices have several interacting state variables: $\dot{x} = f(x)$ with $x$ a vector. For the linear case $\dot{x} = A x$, the eigenvalues of the matrix $A$ are the characteristic roots of §6 — diagonalizing $A$ decouples the system into independent pure cases $y' = \lambda y$. This is the native language of modern control theory.
Laplace transforms and transfer functions. The rule $\mathcal{L}\{y'\} = sY(s) - y(0)$ converts a linear ODE in $t$ into algebra in $s$; the characteristic polynomial reappears as the denominator of $H(s)$, its roots as poles. §7's frequency response is that object evaluated at $s = i\omega$.
Nonlinear dynamics and chaos. Phase portraits generalize Figure 6's sign-reading to systems; from three state variables upward, bounded trajectories can stop repeating while never escaping (the Lorenz system) — the mathematical reason weather forecasts decay after a week even with perfect models.
Partial differential equations. Give the unknown a spatial derivative too — heat flow, waves, fluids — and you have a PDE. The workhorse method of lines discretizes space first, turning the PDE into a large ODE system: back to §8, at scale.
Existence and uniqueness (Picard–Lindelöf, §3) is the theorem working silently behind every solver you run, and finite-time blow-up (the red curve in Figure 6) is why "a solution exists" is a statement worth proving at all.

Summary: one screen to keep

you face	picture	exact tool	if that fails
$y' = k\,y$ (rate ∝ amount)	pure exponential field	$y = y_0 e^{kx}$	thread the field numerically: $y_{n+1} = y_n + h f(x_n, y_n)$ Euler $O(h)$ → RK4 $O(h^4)$ implicit for stiff systems
$y' + p(x)\,y = q(x)$	field with attractor	integrating factor $\mu = e^{\int p}$
$y' = g(x)\,h(y)$ separable	phase line, equilibria	two integrals
$y'' + 2\zeta\omega_n y' + \omega_n^2 y = 0$	mode sum / ringing	characteristic roots
driven oscillator	transfer function peak	$\|H\| = 1/\sqrt{(1-\omega^2/\omega_n^2)^2 + (2\zeta\omega/\omega_n)^2}$

An ODE is a contract between the present and the future: the field says what happens next given what is true now. Read the field, solve the canonical cases exactly — and when the algebra has no answer, walk the field numerically.

Regime	Roots \(r\)	Motion
\(\zeta < 1\) underdamped	\(-\zeta\omega_n \pm i\,\omega_d\), \(\omega_d = \omega_n\sqrt{1-\zeta^2}\)	decaying oscillation \(e^{-\zeta\omega_n t}(A\cos\omega_d t + B\sin\omega_d t)\)
\(\zeta = 1\) critically damped	double root \(-\omega_n\)	\((A + B\,t)\,e^{-\omega_n t}\) — fastest return without overshoot
\(\zeta > 1\) overdamped	two negative reals	slow exponential + fast exponential

you face	picture	exact tool	if that fails
\(y' = k\,y\) (rate ∝ amount)	pure exponential field	\(y = y_0 e^{kx}\)	thread the field numerically: \(y_{n+1} = y_n + h f(x_n, y_n)\) Euler \(O(h)\) → RK4 \(O(h^4)\) implicit for stiff systems
\(y' + p(x)\,y = q(x)\)	field with attractor	integrating factor \(\mu = e^{\int p}\)
\(y' = g(x)\,h(y)\) separable	phase line, equilibria	two integrals
\(y'' + 2\zeta\omega_n y' + \omega_n^2 y = 0\)	mode sum / ringing	characteristic roots
driven oscillator	transfer function peak	\(\|H\| = 1/\sqrt{(1-\omega^2/\omega_n^2)^2 + (2\zeta\omega/\omega_n)^2}\)