This chapter introduces numerical methods for solving ordinary differential equations. Many differential equations that arise in science, engineering, finance, and physics cannot be solved by symbolic techniques. Even when a closed-form solution exists, numerical simulation often provides the fastest and most practical route to understanding how a system behaves. Students will learn to approximate solutions to initial value problems by computing step-by-step trajectories through the solution space.

The methods developed here — Euler's method, the improved Euler method, and the Runge-Kutta family — represent some of the most important and widely used algorithms in all of applied mathematics. Students who have completed the calculus prerequisites for this course and have some exposure to differential equations will find that numerical ODE methods bring together nearly every idea studied so far: function evaluation, approximation, iteration, error analysis, convergence, and stability. The chapter closes with a discussion of systems of differential equations and stiff equations, preparing students for Chapter 13 and the capstone chapter that follows.

Differential equations are the language of change. Every physical process that evolves over time — the temperature of a cooling object, the velocity of a falling body, the voltage across a capacitor, the concentration of a drug in the bloodstream, the spread of an infectious disease, the trajectory of a spacecraft — can be modeled by a differential equation or a system of differential equations.

Symbolic techniques for solving differential equations cover important cases: linear equations with constant coefficients, separable equations, exact equations, and a handful of others. But most differential equations arising from real scientific and engineering problems are nonlinear, coupled, or otherwise resistant to closed-form analysis. Numerical methods make it possible to simulate the behavior of these systems to whatever accuracy a problem demands.

Understanding numerical ODE methods also builds insight into differential equations themselves. Watching an approximation trace out a solution curve step by step gives an intuition for how solutions behave that symbolic formulas sometimes obscure.

A first-order ordinary differential equation expresses a relationship between an unknown function and its derivative: \( \dfrac{dy}{dt} = f(t, y) \). The function \( f(t,y) \) specifies the rate of change of \( y \) at every point in the \( (t,y) \)-plane. Solving the equation means finding a function \( y(t) \) whose derivative satisfies this relationship everywhere on an interval of interest.

When \( f(t,y) \) is simple enough — say, \( f(t,y) = ay \) or \( f(t,y) = g(t) \) — symbolic integration yields an exact formula. The exponential growth equation \( y' = ay \) has the exact solution \( y(t) = y_0 e^{at} \). The equation \( y' = \sin(t) \) can be integrated directly.

But consider \( f(t,y) = \sin(y) + t^2 y \), or a logistic model with a seasonal forcing term, or any coupled system describing the interaction of predator and prey populations. No general technique solves these symbolically. Even when numerical analysts can prove that a solution exists and is unique — the Picard–Lindelöf theorem guarantees this under mild conditions on \( f \) — the proof does not supply a formula.

Numerical methods provide a workaround: instead of finding the exact solution function \( y(t) \), we compute a sequence of approximate values

at a discrete set of time points \( t_0 < t_1 < t_2 < \cdots \). These values trace out an approximate solution curve. If the method is well chosen and the step size is small enough, the approximation can be made as accurate as required.

This philosophy — approximate a continuous problem by a discrete one — runs throughout numerical methods. What is new in this chapter is that the approximation must be built step by step, because each approximate value depends on the previous one.

An initial value problem (IVP) consists of:

1. A first-order differential equation: \( \dfrac{dy}{dt} = f(t,y) \)
2. An initial condition: \( y(t_0) = y_0 \)

The initial condition pins down a specific solution from among all the functions satisfying the differential equation. Without it, the general solution contains an arbitrary constant, and we have infinitely many possibilities.

For numerical work, the initial condition provides the starting point. We know \( y_0 = y(t_0) \) exactly (or as precisely as the problem permits), and we advance forward in time, computing \( y_1, y_2, \ldots \) at successive time values \( t_1, t_2, \ldots \).

The step size \( h \) is the increment between time values:

A uniform step size is simplest and sufficient for most introductory problems. Adaptive methods, which vary \( h \) automatically, are introduced briefly in Section 12.7.

A second-order equation such as \( y'' = f(t, y, y') \) with initial conditions \( y(t_0) = y_0 \) and \( y'(t_0) = v_0 \) can always be rewritten as a system of two first-order equations. This reduction is important: every method in this chapter applies to first-order equations, so higher-order problems and coupled systems are handled by the same framework after this reduction. Section 12.10 treats systems in detail.

Euler's method is the oldest and simplest numerical ODE solver. It is not the most accurate method studied in this chapter, but it is the most transparent: every idea needed to understand more sophisticated methods — discrete steps, local error, global accumulation, step size control — appears in its simplest form here.

Suppose we know \( y(t_n) \) and we want to estimate \( y(t_{n+1}) \) where \( t_{n+1} = t_n + h \). The Taylor expansion of \( y \) around \( t_n \) gives:

Euler's method keeps only the first two terms, replacing the exact derivative \( y'(t_n) \) with the known quantity \( f(t_n, y_n) \):

Geometrically, Euler's method moves along the tangent line to the solution curve at \( (t_n, y_n) \) for a distance \( h \) in the \( t \)-direction. At the new point, it recomputes the slope using \( f(t_{n+1}, y_{n+1}) \) and takes the next step.

Euler's method uses the slope at the left endpoint of each interval to step forward. A natural improvement is to use a better estimate of the slope — specifically, an average of the slope at the beginning and end of the interval.

The difficulty is that we do not know \( y(t_{n+1}) \) yet, so we cannot directly compute the slope at the right endpoint. The improved Euler method (also called Heun's method) resolves this with a two-stage process.

Use Euler's method to produce a preliminary estimate of \( y_{n+1} \):

Average the slopes at the beginning and the predicted end:

This is a predictor-corrector scheme. The predictor gets us a rough endpoint; the corrector uses that endpoint to improve the slope estimate.

The improved Euler method is equivalent to using the trapezoidal rule to approximate the integral in the exact solution formula:

The trapezoidal rule applied to the integrand gives exactly the averaging formula above. Since the trapezoidal rule has local error \( O(h^3) \) compared to the rectangle rule's \( O(h^2) \), the improved Euler method has local truncation error of order \( O(h^3) \) and global error of order \( O(h^2) \). This is the definition of a second-order method: halving the step size reduces the global error by a factor of four.

The Euler and improved Euler methods belong to a broad and important family called Runge-Kutta methods. These methods approximate the solution by taking weighted combinations of function evaluations within each step. No derivatives of \( f \) are required — only evaluations of \( f \) at strategically chosen points.

The key idea is to match terms in a Taylor expansion of the true solution by choosing evaluation points and weights carefully. A method that matches through the term \( h^p \) in the Taylor expansion is called a method of order \( p \):

• Euler's method: order 1 (matches through \( h^1 \))
• Improved Euler: order 2 (matches through \( h^2 \))
• Classical RK4: order 4 (matches through \( h^4 \))

Higher order means that each step is a better approximation of the true solution increment. For smooth equations, higher order dramatically reduces the number of steps needed to achieve a given accuracy.

A general explicit Runge-Kutta method at each step computes several stage values (slopes at intermediate points within the interval \( [t_n, t_{n+1}] \)), then combines them in a weighted sum:

where each \( k_i \) is a slope value computed at some interior point. The specific choices of interior points and weights determine the method's order and accuracy properties.

The classical fourth-order Runge-Kutta method — almost universally referred to as RK4 — is one of the most widely used algorithms in all of numerical computation. It achieves fourth-order accuracy using four function evaluations per step, which is an excellent trade-off between accuracy and cost for smooth problems.

At each step from \( t_n \) to \( t_{n+1} = t_n + h \), compute:

Then advance the solution:

The weights \( 1/6,\, 2/6,\, 2/6,\, 1/6 \) mirror Simpson's rule from numerical integration — and for good reason: RK4 is essentially applying Simpson's rule to approximate the integral of \( f(t, y(t)) \) over each interval.

• \( k_1 \): slope at the left endpoint (same as Euler's method).
• \( k_2 \): slope at the midpoint, estimated using a half-step with slope \( k_1 \).
• \( k_3 \): slope at the midpoint again, re-estimated using slope \( k_2 \) (a refinement).
• \( k_4 \): slope at the right endpoint, estimated using slope \( k_3 \) for a full step.

The method gives extra weight (factor of 2) to the midpoint slopes because midpoint estimates tend to be more accurate than endpoint estimates.

The choice of step size \( h \) is one of the most important decisions in numerical ODE solving. Step size affects accuracy, computational cost, and sometimes stability.

For a method of order \( p \), the local truncation error (LTE) is the error introduced in a single step, assuming the previous value is exact:

for some constant \( C \) that depends on derivatives of \( y \). The local error shrinks rapidly as \( h \) decreases: a method of order 4 has local error proportional to \( h^5 \), so halving \( h \) reduces the local error by a factor of \( 2^5 = 32 \).

The global truncation error (GTE) at a fixed final time \( T = t_0 + Nh \) accumulates over all \( N = (T-t_0)/h \) steps. Since there are \( N \propto 1/h \) steps each contributing \( O(h^{p+1}) \) local error, the global error is:

This is why the order of a method refers to the exponent in the global error. For a method of order \( p \):

A good starting strategy is to run the method with some step size \( h \), then run it again with \( h/2 \), and compare the results. If they agree to the desired number of digits, the step size is likely adequate. If not, halve again and repeat. This is called step doubling or Richardson extrapolation when used systematically.

In practice, the appropriate step size depends on:

• The desired accuracy (how many correct digits are needed).
• The smoothness of \( f(t,y) \) (sharper changes require smaller steps).
• The behavior of the solution (rapidly varying solutions require smaller steps).
• Computational budget (smaller \( h \) means more steps and more time).

As a rule, prefer higher-order methods over tiny step sizes. Using RK4 with \( h = 0.1 \) is generally more efficient than using Euler's method with \( h = 0.001 \) for the same accuracy target.

Understanding the distinction between local and global error is essential for using and comparing ODE solvers reliably.

The local truncation error measures how much the numerical method deviates from the true solution over a single step, assuming the starting value is exact. It is a property of the method and the equation — it does not account for errors that have accumulated from previous steps.

For RK4, the local truncation error is:

for some \( \xi_n \in [t_n, t_{n+1}] \). The fifth derivative of the exact solution appears because RK4 is exact for polynomials of degree up to four.

The global truncation error at time \( T \) is:

where \( y_N \) is the numerical approximation after \( N \) steps. The global error accumulates all local errors, modified by how errors propagate through the equation. If the equation is stable (errors do not grow), the global error is roughly the sum of local errors:

For unstable equations, errors can grow exponentially, and the relationship between local and global error becomes more complicated.

To check that a method is behaving as expected, compute numerical solutions with step sizes \( h \) and \( h/2 \), and compare errors against the exact solution (when known). If the method is order \( p \), the error should decrease by a factor of roughly \( 2^p \) when \( h \) is halved:

This verification is useful for confirming that an implementation is correct and that the step size is in the range where the asymptotic error behavior holds.

Accuracy and stability are related but distinct concerns. A method can be accurate on one equation and completely unreliable on another, even with the same step size.

The standard analysis of numerical stability uses the linear test equation:

where \( \lambda \) is a constant (possibly complex). The exact solution is \( y(t) = e^{\lambda t} \).

If \( \text{Re}(\lambda) < 0 \), the exact solution decays to zero. A numerical method is absolutely stable for a particular value of \( h\lambda \) if the numerical solution also decays to zero (or at least does not grow) for that value.

Applying Euler's method to the test equation gives:

After \( n \) steps:

For the numerical solution to remain bounded, we need \( |1+h\lambda| \leq 1 \). For real \( \lambda < 0 \), this requires:

RK4 has a larger stability region than Euler's method. For the test equation, the RK4 update is:

This is the partial sum of the Taylor series for \( e^{h\lambda} \). For real \( \lambda < 0 \), RK4 is stable for \( h|\lambda| \leq 2.785 \) approximately — a larger interval than Euler's \( h|\lambda| \leq 2 \).

Many real systems involve multiple interacting quantities, each governed by its own differential equation. The equations are coupled: the rate of change of each variable depends on the current values of all variables.

A system of \( m \) first-order equations takes the form:

where \( \mathbf{y} = (y_1, y_2, \ldots, y_m)^T \) is a vector of unknowns, \( \mathbf{y}_0 \) is a vector of initial values, and \( \mathbf{f} = (f_1, f_2, \ldots, f_m)^T \) is a vector-valued function.

Every ODE method in this chapter applies without essential modification: simply replace the scalar \( y \) and scalar \( f \) with the vector \( \mathbf{y} \) and vector \( \mathbf{f} \). The formulas look the same; the computation requires evaluating all components simultaneously.

A second-order equation \( y'' = g(t,y,y') \) can be reduced to a system by introducing a new variable \( v = y' \):

This is a system of two first-order equations in the unknowns \( y \) and \( v \). Initial conditions \( y(t_0) = y_0 \) and \( y'(t_0) = v_0 \) specify the starting point.

Some differential equations pose a special challenge for explicit methods: even though the solution we care about is well-behaved and slowly changing, the equation also admits rapidly decaying transient solutions. These equations are called stiff.

Consider the system \( \dfrac{d\mathbf{y}}{dt} = A\mathbf{y} \), where the matrix \( A \) has eigenvalues with very different magnitudes: one eigenvalue might be \( -1 \) (slow decay) and another \( -1000 \) (very fast decay). The fast-decaying component vanishes quickly and is irrelevant to the long-term solution. But explicit methods must use a step size \( h \leq 2/1000 = 0.002 \) to remain stable, even if all we want to resolve is the slow component at large times.

This mismatch between the time scales of stability and the time scales of interest is the defining feature of stiffness. A stiff equation requires tiny steps for stability even when the solution itself is slowly varying.

Stiff equations arise frequently in:

• Chemical kinetics: Reactions with very different rate constants (some reactions complete in nanoseconds while others take hours).
• Circuit simulation: Electrical circuits with components operating at very different time scales (capacitors and inductors with extreme values).
• Structural dynamics: Mechanical systems with stiff springs or elastic modes at very different frequencies.
• Biological systems: Enzyme kinetics, gene regulation networks, and pharmacokinetic models.

Implicit methods (such as the backward Euler method or implicit Runge-Kutta methods) have much larger stability regions and can take large steps on stiff equations without instability. The trade-off is that implicit methods require solving a system of (possibly nonlinear) equations at each step — more computational work per step, but far fewer steps overall.

Backward Euler method:

Notice that \( y_{n+1} \) appears on both sides. This requires solving for \( y_{n+1} \) at each step, often using Newton's method. For linear equations, this is straightforward. For nonlinear equations, it requires iteration.

The backward Euler method is unconditionally stable for the test equation with \( \text{Re}(\lambda) < 0 \) (regardless of step size), making it suitable for stiff problems where explicit methods would require impractically small steps.

The numerical ODE methods developed in this chapter are not abstract exercises. They are computational engines that simulate the behavior of real systems across every area of science, engineering, medicine, finance, and data science.

Students learning to solve ODEs numerically encounter a predictable set of misconceptions and errors. Recognizing these problems early prevents compounding mistakes across longer computations.

This chapter introduced the fundamental numerical methods for solving ordinary differential equations.