directionField = StreamPlot[{1, x - y}, {x, -0.5, 3}, {y, -2, 4},
StreamStyle -> "PinDart", Frame -> True,
FrameLabel -> {"x", "y"},
PlotLabel -> Style["Direction field of y' = x - y", 13]];
Export["02_ODEs_and_IVPs/direction-field.pdf", directionField];
Export["02_ODEs_and_IVPs/direction-field.svg", directionField];2 Ordinary Differential Equations and Initial Value Problems
2.1 What Is a Differential Equation?
You already know the derivative. Given a function \(y(x)\), its derivative \(y'(x)\) measures how rapidly \(y\) changes as \(x\) increases: it is the slope of the curve at each point. What you have not yet encountered is an equation where the unknown is itself a function, and the equation relates that function to its own derivative.
A differential equation is exactly such an equation. The central class for this book is the first-order ordinary differential equation:
\[\frac{dy}{dx} = f(x, y)\]
Here \(f\) is a given function of two variables — the slope prescription — and \(y(x)\) is the unknown function we seek. The solution is not a number but a curve: one whose slope \(dy/dx\) at every point \(x\) exactly equals \(f(x, y(x))\).
This is a fundamentally different kind of equation from anything in calculus so far. An algebraic equation such as \(x^2 = 4\) asks: for which numbers \(x\) is this statement true? The answer is a finite set. A differential equation asks: for which function \(y(x)\) is this slope constraint true at every \(x\) simultaneously? The answer is a family of curves.
A direction field makes this slope constraint visible. At each point \((x, y)\) in the plane, draw a short arrow with slope \(f(x, y)\). Any solution curve \(y(x)\) must be tangent to the arrow at every point it passes through. The figure below uses \(f(x, y) = x - y\). At any point above the line \(y = x\), the slope \(x - y\) is negative and the curve must descend. Below \(y = x\), the slope is positive and the curve must rise. All solutions are drawn toward the long-run trend \(y \approx x - 1\).
Trace any arrow in the figure. A solution starting far above the diagonal curves steeply downward. As \(y\) approaches \(x - 1\), the slope approaches zero and the curve levels off. Every solution, regardless of starting point, exhibits this convergence — a qualitative fact you can read directly from the direction field, without solving a single equation.
The direction field is the geometric face of the differential equation. Its algebraic face is the formula \(f(x, y)\). The function \(y(x)\) we seek must make these two faces agree everywhere.
2.2 Initial Value Problems: Pinning Down a Unique Solution
A single differential equation typically has infinitely many solutions. Consider \(y' = y\): we require that the slope of \(y\) at every \(x\) equals its current value there. It is easy to verify that for any constant \(C\), the function \(y(x) = C\, e^x\) satisfies this. Differentiating: the derivative of \(C\, e^x\) is \(C\, e^x\), which equals \(y(x)\). The equation holds for every \(C\).
\[y(x) = C\, e^x, \quad C \in \mathbb{R}\]
The code below plots this family for several values of \(C\). Every curve is a valid solution to \(y' = y\). They fill the plane without crossing.
familyPlot = Plot[Table[c * Exp[x], {c, {-2, -1, -0.5, 0.5, 1, 2, 3}}],
{x, -0.5, 2},
PlotRange -> {-5, 22}, Frame -> True,
FrameLabel -> {"x", "y"},
PlotLabel -> Style["Family of solutions: y = C Exp[x]", 13]];
Export["02_ODEs_and_IVPs/family.pdf", familyPlot];
Export["02_ODEs_and_IVPs/family.svg", familyPlot];
To select one solution from this family, we supply an additional piece of information: the value of \(y\) at a specific \(x\). This is the initial condition, and the combination of a differential equation with an initial condition is called an initial value problem, or IVP:
\[\frac{dy}{dx} = f(x, y), \quad y(x_0) = y_0\]
For \(y' = y\) with \(y(0) = 1\): the solution must pass through \((0, 1)\). Since \(y(0) = C\, e^0 = C\), we need \(C = 1\), and the unique solution is \(y(x) = e^x\). The initial condition collapses the infinite family to a single curve.
This uniqueness is guaranteed by the Picard–Lindelöf theorem: if \(f(x, y)\) and the partial derivative \(\partial f / \partial y\) are both continuous near \((x_0, y_0)\), then the IVP has exactly one solution in a neighborhood of \(x_0\). The continuity of \(\partial f / \partial y\) ensures the slope field cannot change so abruptly in the \(y\)-direction that two distinct solution curves might pass through the same point.
The figure below shows four different initial conditions \((x_0, y_0)\) overlaid on the slope field of \(y' = y\). The red dot marks each initial condition; the blue curve is the unique solution it determines. Every starting point selects exactly one curve from the family \(C\, e^x\).
ivpPanels = Table[
With[{x0 = ic[[1]], y0 = ic[[2]]},
Show[
StreamPlot[{1, y}, {x, -0.5, 2.5}, {y, -3, 9},
StreamStyle -> {Gray, Thin}, Frame -> True,
FrameLabel -> {"x", "y"}],
Plot[y0 * Exp[x - x0], {x, -0.5, 2.5},
PlotStyle -> {Thick, Blue}, PlotRange -> {-3, 9}],
Graphics[{Red, PointSize[0.025], Point[{x0, y0}]}],
PlotLabel -> Style[
"x0 = " <> ToString[x0] <> ", y0 = " <> ToString[y0], 11],
ImageSize -> 280
]
],
{ic, {{0, 1}, {0.5, 0.5}, {1, 2}, {0, -1}}}
];
ivpGrid = GraphicsGrid[Partition[ivpPanels, 2], Spacings -> Scaled[0.04]];
Export["02_ODEs_and_IVPs/ivp-panels.pdf", ivpGrid];
Export["02_ODEs_and_IVPs/ivp-panels.svg", ivpGrid];
The direction field shows all solutions simultaneously; the initial condition picks one.
2.3 When Analytical Solutions Exist — and When They Don’t
The Picard–Lindelöf theorem guarantees a solution exists. It says nothing about whether that solution can be written in terms of familiar functions — polynomials, exponentials, logarithms, trigonometric functions, and their finite combinations. For most differential equations arising in practice, no such formula exists.
Two classes with systematic closed-form solutions are worth knowing, because they illuminate what “solvable” means and what it requires.
Separable equations take the form \(y' = g(x) h(y)\). The right side factors into a part depending only on \(x\) and a part depending only on \(y\). Dividing both sides by \(h(y)\) and integrating: the left integral involves only \(y\), the right only \(x\). This separates the equation into two standard integrals that can be solved independently.
Linear first-order equations take the form \(y' + p(x) y = q(x)\). Multiplying through by the integrating factor \(\mu(x) = \exp\!\left(\int p(x)\, dx\right)\) converts the left side into the exact derivative \((\mu y)'\), which can then be integrated directly.
The equation \(y' = y\) is both separable (write \(dy/y = dx\) and integrate both sides) and linear (\(p(x) = -1\), \(q(x) = 0\)), so DSolve[] has no difficulty:
DSolve[y'[x] == y[x], y[x], x]
The output is {y[x] -> C[1] E^x}, confirming the one-parameter family. Adding the initial condition selects the unique solution:
DSolve[{y'[x] == y[x], y[0] == 1}, y[x], x]
Now consider \(y' = \sin(x y)\). The right side mixes \(x\) and \(y\) in a way that is neither separable nor linear — there is no algebraic way to isolate \(x\) from \(y\) or to identify an integrating factor.
DSolve[y'[x] == Sin[x * y[x]], y[x], x]
Wolfram Language returns the equation unevaluated. This is not a software limitation; no computer algebra system can produce a closed form because one does not exist. Yet the solution \(y(x)\) is perfectly well-defined for any initial condition satisfying the Picard–Lindelöf conditions. We can compute it numerically to any desired precision. We simply cannot write it as a formula.
This situation is the rule, not the exception. The equations governing a pendulum with air resistance, a driven electrical circuit with a nonlinear component, or a population with logistic coupling all fall outside the two solvable classes. Numerical methods are the standard tool.
2.4 The Numerical Strategy: Stepping Forward in Time
When a closed form is unavailable, we compute approximate values of \(y\) at a sequence of \(x\)-positions. The strategy is direct.
Suppose we know \(y_n = y(x_n)\) exactly. The differential equation tells us the slope at that point: it is \(f(x_n, y_n)\). If we step a small distance \(h\) to the right to reach \(x_n + h\), the first-order Taylor approximation gives:
\[y_{n+1} \approx y_n + h \cdot f(x_n, y_n)\]
Applying this repeatedly, starting from \((x_0, y_0)\), generates approximate values at \(x_1 = x_0 + h\), \(x_2 = x_0 + 2h\), and so on. Each step uses only the current position and the slope at that position. The approximation assumes the slope is constant over the step interval, which it is not — the slope changes as \(y\) changes. The error from this assumption is proportional to \(h^2\) per step.
The figure below applies this stepping formula to \(y' = y\) with \(y(0) = 1\), whose exact solution is \(e^x\). The blue curve is the exact solution; the red polyline connects the computed points \((x_n, y_n)\). Four step sizes are shown side by side.
hValues = {1.0, 0.5, 0.25, 0.1};
eulerPanels = Table[
Module[{pts},
pts = NestList[{#[[1]] + h, #[[2]] + h * #[[2]]} &, {0., 1.}, Round[2./h]];
Show[
Plot[Exp[x], {x, 0, 2}, PlotStyle -> {Blue, Thick}, PlotRange -> {0, 8}],
ListLinePlot[pts, PlotStyle -> Red, Mesh -> All,
MeshStyle -> {Red, PointSize[0.018]}, PlotRange -> {0, 8}],
Frame -> True, FrameLabel -> {"x", "y"},
PlotLabel -> Style[
"h = " <> ToString[h] <> ", steps = " <> ToString[Round[2./h]], 11],
ImageSize -> 280
]
],
{h, hValues}
];
eulerGrid = GraphicsGrid[Partition[eulerPanels, 2], Spacings -> Scaled[0.04]];
Export["02_ODEs_and_IVPs/euler-grid.pdf", eulerGrid];
Export["02_ODEs_and_IVPs/euler-grid.svg", eulerGrid];
With \(h = 1\), one slope evaluation per unit interval, the approximation at \(x = 2\) overshoots the true value \(e^2 \approx 7.39\) noticeably. With \(h = 0.1\), it stays close. Halving \(h\) reduces the per-step error by a factor of roughly \(4\) (since it is proportional to \(h^2\)), at the cost of twice as many function evaluations. This accuracy-vs-cost tradeoff is central to all numerical ODE methods.
This formula — use the slope at the left endpoint of each interval — is Euler’s method. Chapter 4 derives it formally from the Taylor series and quantifies its error precisely. Every subsequent chapter asks the same question: can we evaluate the slope function at more than one point per step, combine those evaluations cleverly, and achieve much higher accuracy for the same \(h\)?
2.5 Setting Up Our Test Problem for the Rest of the Book
Every method in this book — Euler, RK2, RK3, and RK4 — will be benchmarked against the same initial value problem:
\[\frac{dy}{dx} = y, \quad y(0) = 1\]
The exact solution is \(y(x) = e^x\). This benchmark earns its place for three reasons. First, the exact answer is known, so we can compute the numerical error precisely at every step rather than estimating it. Second, \(e^x\) is smooth on any finite interval and satisfies the Picard–Lindelöf conditions without subtlety, so we never have to worry about existence or uniqueness. Third, the slope function \(f(x, y) = y\) is the simplest possible — just the current value of \(y\) — yet the exponential growth it drives is nontrivial enough to reveal real differences between methods.
We define the test problem in Wolfram Language now. Every subsequent chapter will reload these definitions at the top of its notebook.
f[x_, y_] := yexactSolution[x_] := Exp[x]x0 = 0; y0 = 1; xEnd = 2;Two quick consistency checks: the exact solution at \(x_0\) must equal \(y_0\), and the slope function at the starting point must return the value \(f(0, 1) = 1\) (since \((e^x)'\big|_{x=0} = 1\)).
exactSolution[x0] == y0
f[x0, y0]1
Both return 1. The slope at the starting point equals the function value, which is precisely the self-referential property that drives exponential growth.
Here is the exact solution plotted over the integration interval \([0, 2]\). This is the target curve that every subsequent chapter’s numerical method will try to match.
exactPlot = Plot[exactSolution[x], {x, 0, 2}, Frame -> True,
FrameLabel -> {"x", "y"}, PlotStyle -> {Blue, Thick},
PlotLabel -> Style["Exact solution: y(x) = Exp[x]", 13],
PlotRangePadding -> 0.3];
Export["02_ODEs_and_IVPs/exact-solution.pdf", exactPlot];
Export["02_ODEs_and_IVPs/exact-solution.svg", exactPlot];