ODE Solvers

Why ODE Solvers Matters

A cooling model says how quickly temperature changes at the current temperature, but it does not directly give the temperature an hour later. A differential-equation solver starts from a known state and repeatedly uses those rates of change to approximate the future state. The same idea drives motion simulation, chemical kinetics, and dynamical-system models.

This lesson solves a real initial-value problem, $\mathbf y'(t)=\mathbf f(t,\mathbf y)$ with $\mathbf y(t_0)=\mathbf y_0$. It develops forward Euler, two-stage Heun, classical RK4, and an adaptive Runge–Kutta method using the Bogacki–Shampine $3(2)$ pair. You will trace every stage on a small equation, implement vector states in three languages, and distinguish local error estimates, actual solution error, and numerical stability.

The implementations are explicit methods for nonstiff problems. A right-hand-side function supplies rates for all state components; no analytic solution is required. An analytic solution is used in the sample only to check the computed error independently.

Core Concepts and Terminology

For a scalar example, take $y'=-y$ and $y(0)=1$, whose exact solution is $y(t)=e^{-t}$. At the initial point the slope is $-1$. An Euler step of length $h=1/2$ follows that starting tangent and produces $1-(1/2)=1/2$; the exact value at the new time is $e^{-1/2}\approx0.606531$. Sampling additional slopes can improve the step.

For a system, $\mathbf y$ is a vector of $d$ components and $\mathbf f$ returns a vector of the same length. For example, the oscillator $x''=-x$ becomes $\mathbf y=[x,v]^T$ and $\mathbf f=[v,-x]^T$. All stage-state updates apply componentwise. Explicit methods evaluate known trial states rather than solving a matrix system; they have no pivoting step.

Assume the right-hand side defines a sufficiently smooth solution on the interval of interest. Local Lipschitz continuity in the state gives the usual local uniqueness setting; it does not guarantee a solution remains finite for all time. A solver's finite-value checks cannot establish existence, smoothness, or a model's physical validity.

How ODE Solvers Works, Step by Step

1. Validate the interval, state, and method parameters. Require finite endpoints with a finite representable difference, a nonempty finite initial vector, and a fixed right-hand-side dimension. Fixed solvers take a positive step count. The adaptive solver takes a positive initial step magnitude, positive absolute tolerance, nonnegative relative tolerance, and a nonnegative attempt budget. A validated zero-length interval returns the initial state without evaluating the right-hand side.

2. Choose a step that stays inside the requested interval. The fixed code schedules $N$ steps from the original interval and explicitly sets its final endpoint to $t_1$. The adaptive code chooses the direction from $t_1-t_0$ and clips each trial to the remaining interval. Before any evaluations, it checks that the endpoint advances and that required interior stage times are representable. Failure to form the stage times returns stalled at the last committed state.

3. Construct stage states and evaluate slopes. Here $\mathbf k_i$ always means a derivative, without an included factor of $h$.

Forward Euler uses one slope:

$$ \mathbf k_1=\mathbf f(t,\mathbf y),\qquad \mathbf y_{\mathrm{new}}=\mathbf y+h\mathbf k_1. $$

Heun's method here means the explicit trapezoidal method: predict with Euler, then average the starting and predicted endpoint slopes:

$$ \mathbf k_1=\mathbf f(t,\mathbf y),\quad \mathbf k_2=\mathbf f(t+h,\mathbf y+h\mathbf k_1),\quad \mathbf y_{\mathrm{new}}=\mathbf y+\frac h2(\mathbf k_1+\mathbf k_2). $$

This is an explicit predictor-corrector formula, not the implicit trapezoidal rule. Classical RK4 uses four slopes:

$$ \begin{aligned} \mathbf k_1&=\mathbf f(t,\mathbf y),\\ \mathbf k_2&=\mathbf f(t+h/2,\mathbf y+h\mathbf k_1/2),\\ \mathbf k_3&=\mathbf f(t+h/2,\mathbf y+h\mathbf k_2/2),\\ \mathbf k_4&=\mathbf f(t+h,\mathbf y+h\mathbf k_3),\\ \mathbf y_{\mathrm{new}}&=\mathbf y+\frac h6(\mathbf k_1+2\mathbf k_2+2\mathbf k_3+\mathbf k_4). \end{aligned} $$

Stages two and three occur at the same time but generally use different states. Overwriting the committed state between stages would change the method. The code stores stage vectors separately and constructs each stage from the original state for that trial.

4. For adaptivity, form two estimates using the Bogacki–Shampine pair. With four shared evaluations, define

$$ \begin{aligned} \mathbf k_1&=\mathbf f(t,\mathbf y),\\ \mathbf k_2&=\mathbf f(t+h/2,\mathbf y+h\mathbf k_1/2),\\ \mathbf k_3&=\mathbf f(t+3h/4,\mathbf y+3h\mathbf k_2/4),\\ \mathbf y^{[3]}&=\mathbf y+h(2\mathbf k_1/9+\mathbf k_2/3+4\mathbf k_3/9),\\ \mathbf k_4&=\mathbf f(t+h,\mathbf y^{[3]}),\\ \mathbf y^{[2]}&=\mathbf y+h(7\mathbf k_1/24+\mathbf k_2/4+\mathbf k_3/3+\mathbf k_4/8). \end{aligned} $$

The bracketed superscripts denote orders, not powers. If accepted, the code advances with the third-order estimate. The difference between the two formulas is generally $O(h^3)$ and estimates the lower-order formula's local error; it is not a certified bound on the retained third-order state or the final global error.

5. Normalize the error across components and accept or reject. For each component, use

$$ s_j=\text{atol}+\text{rtol}\max(\lvert y_j\rvert,\lvert y_j^{[3]}\rvert),\qquad E=\max_j\frac{\lvert y_j^{[3]}-y_j^{[2]}\rvert}{s_j}. $$

The programs use scalar tolerances and a maximum norm. Positive atol keeps every scale positive near zero; component-specific tolerances are a useful extension when state units differ. If $E\le1$, commit the time and third-order state and append one trajectory row. If $E\gt1$, reject the trial and leave time, state, and stored history unchanged. Either outcome counts the trial's four evaluations.

6. Adjust the next step magnitude and finish with a reason. For nonzero $E$, multiply the last trial's step magnitude by $\operatorname{clip}(0.9E^{-1/3},0.2,5)$. The cube-root exponent comes from the $O(h^3)$ difference. For $E=0$, use factor $5$; after a rejection, also cap the factor at $1$. Clip the result to the interval still to integrate. An overflowing normalized error is treated as an infinite error and rejected; nonfinite function values, stage states, or checked tolerance scales raise an exception.

flowchart TD
    A["Validate interval, initial vector, and controls"] --> B{"Already at target?"}
    B -->|Yes| R["Return done with committed trajectory"]
    B -->|No| C{"Attempt budget remains?"}
    C -->|No| L["Return limit with latest committed state"]
    C -->|Yes| D["Clip trial endpoint to remaining interval"]
    D --> E{"Stage times representable?"}
    E -->|No| S["Return stalled"]
    E -->|Yes| F["Evaluate four shared stages and both trial states"]
    F --> G["Compute scaled maximum error E"]
    G --> H{"E at most 1?"}
    H -->|Yes| I["Commit time and third-order state; append history"]
    H -->|No| J["Count rejection; retain time, state, and history"]
    I --> K{"Target reached?"}
    K -->|Yes| R
    K -->|No| M["Adjust step magnitude"]
    J --> M
    M --> C

The diagram shows the adaptive driver. Fixed-step methods commit every finite, representable trial and have no embedded error acceptance test. done means the target time was reached, limit means the adaptive attempt budget was exhausted, and stalled means the time grid could no longer support the required stages. None of these statuses by itself certifies a global error bound.

Worked Example with a Full Trace

Solve $y'=-y$, $y(0)=1$ on $[0,1]$. For the readable hand trace, use $h=1/2$, giving two steps for each fixed method. The exact endpoint is $e^{-1}\approx0.3678794412$.

Euler maps a starting value $q$ to $q/2$. Heun's two slopes are $-q$ and $-q/2$, so it maps $q$ to $5q/8$. Both complete trajectories are shown in the final table below.

RK4's complete stage trace is:

Each RK4 step multiplies its current value by $233/384$, so the final state is $(233/384)^2\approx0.3681708442$. Fractions describe exact-arithmetic stage values; the programs store floating-point approximations.

For a complete short adaptive trace, use initial step magnitude $1/2$, $\text{atol}=10^{-2}$, $\text{rtol}=0$, and attempt budget $10$. Its stages on the first trial are $k_1=-1$, $k_2=-3/4$, $k_3=-23/32$, and $k_4=-29/48$. They produce $y^{[3]}=29/48$ and $y^{[2]}=463/768$.

The second trial's slopes are the first trial's slopes multiplied by $29/48$: $-29/48$, $-29/64$, $-667/1536$, and $-841/2304$. After the first acceptance the controller proposes growth, but clipping leaves exactly $1/2$ to reach the target. The run ends with two accepted steps, no rejections, and eight evaluations.

All four hand trajectories, with every committed state, are:

For a rejection contrast, change only the adaptive absolute tolerance to $10^{-3}$. The first trial now has $E=125/96\approx1.3020833$ and is rejected. With an attempt budget of one, the returned status is limit, time remains $0$, state remains $1$, history contains only the initial row, and the evaluation count is $4$.

Implementation in Python, C++, and Java

Each panel is a complete program using ordinary state arrays, nested arrays for Runge–Kutta coefficients and stages, and trajectory rows of the form $[t,y_0,\ldots,y_{d-1}]$. A single tableau-driven step routine implements the four formulas. These explicit methods require no linear-system solve or matrix pivoting; implicit methods discussed later would require a different step routine.

fixed accepts a right-hand-side callback, start and target times, initial vector, positive step count, and method string euler, heun, or rk4. adaptive uses the Bogacki–Shampine pair with the tolerance and attempt controls described above. Backward integration is supported by reversing the time direction; h0 remains a positive magnitude. The callback must return a finite vector of the original dimension and should represent a deterministic rate law. It receives a separate stage-state copy.

The result stores final time, state, full accepted trajectory, accepted and rejected counts, evaluation count, and status. Rejected candidates are not stored. These implementations intentionally do not reuse the last stage as the next first stage, so each adaptive attempt costs exactly four right-hand-side evaluations. A production Bogacki–Shampine implementation can exploit that reuse.

The executable sample uses the same decay problem on $[0,1]$, but four fixed steps rather than the hand trace's two. The adaptive sample uses initial magnitude $0.5$, absolute tolerance $10^{-7}$, relative tolerance $10^{-6}$, and attempt budget $10000$. Its independent check compares with $e^{-1}$; local tolerances are not substituted for actual endpoint error.

import math
from dataclasses import dataclass


def finite(v):
    if not math.isfinite(v):
        raise ValueError("nonfinite input, stage, or arithmetic result")
    return v


def validate(t0, t1, y0):
    finite(t0)
    finite(t1)
    finite(t1-t0)
    if not y0:
        raise ValueError("nonempty initial state required")
    return [finite(v) for v in y0]


def tableau(method):
    if method == "euler":
        return [0.0], [[0.0]], [1.0], [1.0]
    if method == "heun":
        return [0.0,1.0], [[0.0,0.0],[1.0,0.0]], [0.5,0.5], [0.5,0.5]
    if method == "rk4":
        return ([0.0,0.5,0.5,1.0],
                [[0.,0.,0.,0.],[0.5,0.,0.,0.],[0.,0.5,0.,0.],[0.,0.,1.,0.]],
                [1/6,1/3,1/3,1/6], [1/6,1/3,1/3,1/6])
    if method == "bs23":
        return ([0.0,0.5,0.75,1.0],
                [[0.,0.,0.,0.],[0.5,0.,0.,0.],[0.,0.75,0.,0.],[2/9,1/3,4/9,0.]],
                [2/9,1/3,4/9,0.], [7/24,1/4,1/3,1/8])
    raise ValueError("unknown method")


def times_ok(t, end, method):
    if not math.isfinite(end) or end == t:
        return False
    h = end-t
    lo, hi = min(t,end), max(t,end)
    mid = t+0.5*h
    if method in ("rk4","bs23") and not lo < mid < hi:
        return False
    if method == "bs23":
        q = t+0.75*h
        if not min(mid,end) < q < max(mid,end):
            return False
    return True


def rk_step(f, t, y, end, method):
    c,a,b,low = tableau(method)
    h, stages = end-t, []
    for k in range(len(c)):
        state = []
        for j in range(len(y)):
            total = 0.0
            for i in range(k):
                total += a[k][i]*stages[i][j]
            state.append(finite(y[j]+h*total))
        time = end if c[k] == 1.0 else finite(t+c[k]*h)
        slope = f(time,state.copy())
        if len(slope) != len(y):
            raise ValueError("right-hand-side dimension mismatch")
        stages.append([finite(v) for v in slope])
    high, lower = [], []
    for j in range(len(y)):
        s1, s2 = 0.0, 0.0
        for k in range(len(c)):
            s1 += b[k]*stages[k][j]
            s2 += low[k]*stages[k][j]
        high.append(finite(y[j]+h*s1))
        lower.append(finite(y[j]+h*s2))
    return high,lower,len(c)


@dataclass
class Result:
    t: float
    y: list
    path: list
    accepted: int
    rejected: int
    evaluations: int
    status: str


def fixed(f, t0, t1, y0, n, method):
    y = validate(t0,t1,y0)
    if method not in ("euler","heun","rk4") or not isinstance(n,int) or n < 1:
        raise ValueError("valid fixed method and positive step count required")
    t, path, calls = t0, [[t0]+y.copy()], 0
    if t0 == t1:
        return Result(t,y,path,0,0,0,"done")
    for i in range(n):
        end = t1 if i == n-1 else finite(t0+(t1-t0)*((i+1)/n))
        if not times_ok(t,end,method):
            return Result(t,y,path,i,0,calls,"stalled")
        y,_,used = rk_step(f,t,y,end,method)
        calls += used
        t = end
        path.append([t]+y.copy())
    return Result(t,y,path,n,0,calls,"done")


def adaptive(f, t0, t1, y0, h0=0.5, atol=1e-7, rtol=1e-6, limit=10000):
    y = validate(t0,t1,y0)
    if (not math.isfinite(h0) or h0 <= 0.0 or not math.isfinite(atol) or atol <= 0.0
            or not math.isfinite(rtol) or rtol < 0.0 or not isinstance(limit,int) or limit < 0):
        raise ValueError("invalid initial step, tolerances, or attempt budget")
    t, path, accepted, rejected, calls = t0, [[t0]+y.copy()], 0, 0, 0
    if t0 == t1:
        return Result(t,y,path,0,0,0,"done")
    direction, size = (1.0 if t1 > t0 else -1.0), h0
    for _ in range(limit):
        remaining = abs(t1-t)
        size = min(size,remaining)
        end = t1 if size == remaining else finite(t+direction*size)
        if not times_ok(t,end,"bs23"):
            return Result(t,y,path,accepted,rejected,calls,"stalled")
        high,low,used = rk_step(f,t,y,end,"bs23")
        calls += used
        error = 0.0
        for old,new,other in zip(y,high,low):
            scale = finite(atol+rtol*max(abs(old),abs(new)))
            error = max(error,abs(finite(new-other))/scale)
        h = abs(end-t)
        if error <= 1.0:
            t, y = end, high
            accepted += 1
            path.append([t]+y.copy())
            if t == t1:
                return Result(t,y,path,accepted,rejected,calls,"done")
        else:
            rejected += 1
        factor = 5.0 if error == 0.0 else max(0.2,min(5.0,0.9*error**(-1.0/3.0)))
        if error > 1.0:
            factor = min(1.0,factor)
        size = min(h*factor,abs(t1-t))
    return Result(t,y,path,accepted,rejected,calls,"limit")


def main():
    def f(t,y):
        return [-y[0]]

    exact = math.exp(-1.0)
    results = [fixed(f,0.,1.,[1.],4,method) for method in ["euler","heun","rk4"]]
    results.append(adaptive(f,0.,1.,[1.]))
    names, bounds = ["euler","heun","rk4","adaptive-bs23"], [0.06,0.005,2e-5,2e-6]
    ok = True
    for name,r,bound in zip(names,results,bounds):
        error = abs(r.y[0]-exact)
        ok = ok and r.status == "done" and r.t == 1.0 and error <= bound
        print(f"{name}: y(1) = {r.y[0]:.6f}, actual error = {error:.6f}")
    r = results[-1]
    print(f"adaptive status = {r.status}")
    print(f"adaptive accepted = {r.accepted}, rejected = {r.rejected}, evaluations = {r.evaluations}")
    print(f"check = {'PASS' if ok else 'FAIL'}")


if __name__ == "__main__":
    main()

#include <algorithm>
#include <cmath>
#include <functional>
#include <iomanip>
#include <iostream>
#include <stdexcept>
#include <string>
#include <vector>
using Vec=std::vector<double>;using Matrix=std::vector<Vec>;
using Fn=std::function<Vec(double,Vec)>;
double finite_value(double v) {
    if(!std::isfinite(v)) throw std::invalid_argument("nonfinite input, stage, or arithmetic result");
    return v;
}
Vec validate(double t0,double t1,const Vec& y0) {
    finite_value(t0);finite_value(t1);finite_value(t1-t0);
    if(y0.empty()) throw std::invalid_argument("nonempty initial state required");
    for(double v:y0)finite_value(v);
    return y0;
}
struct Tableau {Vec c;Matrix a;Vec b,low;};
Tableau tableau(const std::string& method) {
    if(method=="euler")return { {0.}, { {0.} }, {1.}, {1.} };
    if(method=="heun")return { {0.,1.}, { {0.,0.}, {1.,0.} }, {.5,.5}, {.5,.5} };
    if(method=="rk4")return { {0.,.5,.5,1.},
        { {0.,0.,0.,0.}, {.5,0.,0.,0.}, {0.,.5,0.,0.}, {0.,0.,1.,0.} },
        {1./6,1./3,1./3,1./6}, {1./6,1./3,1./3,1./6} };
    if(method=="bs23")return { {0.,.5,.75,1.},
        { {0.,0.,0.,0.}, {.5,0.,0.,0.}, {0.,.75,0.,0.}, {2./9,1./3,4./9,0.} },
        {2./9,1./3,4./9,0.}, {7./24,1./4,1./3,1./8} };
    throw std::invalid_argument("unknown method");
}
bool times_ok(double t,double end,const std::string& method) {
    if(!std::isfinite(end) || end==t)return false;
    double h=end-t,lo=std::min(t,end),hi=std::max(t,end),mid=t+.5*h;
    if((method=="rk4" || method=="bs23") && !(lo<mid && mid<hi))return false;
    if(method=="bs23") {
        double q=t+.75*h;if(!(std::min(mid,end)<q && q<std::max(mid,end)))return false;
    }
    return true;
}
struct Trial {Vec high,low;int evaluations;};
Trial rk_step(const Fn& f,double t,const Vec& y,double end,const std::string& method) {
    auto tab=tableau(method);double h=end-t;Matrix stages;int d=static_cast<int>(y.size());
    int s=static_cast<int>(tab.c.size());
    for(int k=0;k<s;++k) {
        Vec state(d);
        for(int j=0;j<d;++j) {
            double total=0;for(int i=0;i<k;++i)total+=tab.a[k][i]*stages[i][j];
            state[j]=finite_value(y[j]+h*total);
        }
        double time=tab.c[k]==1. ? end : finite_value(t+tab.c[k]*h);Vec slope=f(time,state);
        if(slope.size()!=y.size())throw std::invalid_argument("right-hand-side dimension mismatch");
        for(double v:slope)finite_value(v);
        stages.push_back(slope);
    }
    Vec high(d),low(d);
    for(int j=0;j<d;++j) {
        double s1=0,s2=0;for(int k=0;k<s;++k) {s1+=tab.b[k]*stages[k][j];s2+=tab.low[k]*stages[k][j];}
        high[j]=finite_value(y[j]+h*s1);low[j]=finite_value(y[j]+h*s2);
    }
    return {high,low,s};
}
struct Result {double t;Vec y;Matrix path;int accepted,rejected;long long evaluations;std::string status;};
Result fixed(const Fn& f,double t0,double t1,const Vec& y0,int n,const std::string& method) {
    Vec y=validate(t0,t1,y0);
    if((method!="euler" && method!="heun" && method!="rk4") || n<1)
        throw std::invalid_argument("valid fixed method and positive step count required");
    double t=t0;Matrix path;Vec row={t};row.insert(row.end(),y.begin(),y.end());path.push_back(row);
    long long calls=0;if(t0==t1)return {t,y,path,0,0,0,"done"};
    for(int i=0;i<n;++i) {
        double end=i==n-1 ? t1 : finite_value(t0+(t1-t0)*(static_cast<double>(i+1)/n));
        if(!times_ok(t,end,method))return {t,y,path,i,0,calls,"stalled"};
        auto trial=rk_step(f,t,y,end,method);y=trial.high;calls+=trial.evaluations;t=end;
        row={t};row.insert(row.end(),y.begin(),y.end());path.push_back(row);
    }
    return {t,y,path,n,0,calls,"done"};
}
Result adaptive(const Fn& f,double t0,double t1,const Vec& y0,double h0=.5,double atol=1e-7,double rtol=1e-6,int limit=10000) {
    Vec y=validate(t0,t1,y0);
    if(!std::isfinite(h0) || h0<=0 || !std::isfinite(atol) || atol<=0 || !std::isfinite(rtol) || rtol<0 || limit<0)
        throw std::invalid_argument("invalid initial step, tolerances, or attempt budget");
    double t=t0;Vec row={t};row.insert(row.end(),y.begin(),y.end());Matrix path={row};
    int accepted=0,rejected=0;long long calls=0;
    if(t0==t1)return {t,y,path,0,0,0,"done"};
    double direction=t1>t0 ? 1. : -1.,size=h0;
    for(int attempt=0;attempt<limit;++attempt) {
        double remaining=std::abs(t1-t);size=std::min(size,remaining);
        double end=size==remaining ? t1 : finite_value(t+direction*size);
        if(!times_ok(t,end,"bs23"))return {t,y,path,accepted,rejected,calls,"stalled"};
        auto trial=rk_step(f,t,y,end,"bs23");calls+=trial.evaluations;double error=0;
        for(std::size_t j=0;j<y.size();++j) {
            double scale=finite_value(atol+rtol*std::max(std::abs(y[j]),std::abs(trial.high[j])));
            error=std::max(error,std::abs(finite_value(trial.high[j]-trial.low[j]))/scale);
        }
        double h=std::abs(end-t);
        if(error<=1.) {
            t=end;y=trial.high;++accepted;row={t};row.insert(row.end(),y.begin(),y.end());path.push_back(row);
            if(t==t1)return {t,y,path,accepted,rejected,calls,"done"};
        } else ++rejected;
        double factor=error==0 ? 5. : std::max(.2,std::min(5.,.9*std::pow(error,-1./3)));
        if(error>1.)factor=std::min(1.,factor);
        size=std::min(h*factor,std::abs(t1-t));
    }
    return {t,y,path,accepted,rejected,calls,"limit"};
}
int main() {
    Fn f=[](double,Vec y){return Vec{-y[0]};};double exact=std::exp(-1.);
    std::vector<Result> results;
    for(const auto& method:std::vector<std::string>{"euler","heun","rk4"})results.push_back(fixed(f,0,1,{1.},4,method));
    results.push_back(adaptive(f,0,1,{1.}));
    std::vector<std::string> names={"euler","heun","rk4","adaptive-bs23"};Vec bounds={.06,.005,2e-5,2e-6};bool ok=true;
    std::cout<<std::fixed<<std::setprecision(6);
    for(int i=0;i<4;++i) {
        const auto& r=results[i];double error=std::abs(r.y[0]-exact);
        ok=ok && r.status=="done" && r.t==1. && error<=bounds[i];
        std::cout<<names[i]<<": y(1) = "<<r.y[0]<<", actual error = "<<error<<'\n';
    }
    const auto& r=results.back();std::cout<<"adaptive status = "<<r.status<<'\n';
    std::cout<<"adaptive accepted = "<<r.accepted<<", rejected = "<<r.rejected<<", evaluations = "<<r.evaluations<<'\n';
    std::cout<<"check = "<<(ok?"PASS":"FAIL")<<'\n';return 0;
}

import java.util.ArrayList;
import java.util.Locale;
public class Main {
    interface Fn {double[] apply(double t,double[] y);}
    static double finite(double v) {
        if(!Double.isFinite(v))throw new IllegalArgumentException("nonfinite input, stage, or arithmetic result");
        return v;
    }
    static double[] validate(double t0,double t1,double[] y0) {
        finite(t0);finite(t1);finite(t1-t0);
        if(y0.length==0)throw new IllegalArgumentException("nonempty initial state required");
        for(double v:y0)finite(v);return y0.clone();
    }
    static class Tableau {
        double[] c,b,low;double[][] a;
        Tableau(double[] c,double[][] a,double[] b,double[] low){this.c=c;this.a=a;this.b=b;this.low=low;}
    }
    static Tableau tableau(String method) {
        if(method.equals("euler"))return new Tableau(new double[]{0.},new double[][]{ {0.} },new double[]{1.},new double[]{1.});
        if(method.equals("heun"))return new Tableau(new double[]{0.,1.},new double[][]{ {0.,0.}, {1.,0.} },new double[]{.5,.5},new double[]{.5,.5});
        if(method.equals("rk4"))return new Tableau(new double[]{0.,.5,.5,1.},
            new double[][]{ {0.,0.,0.,0.}, {.5,0.,0.,0.}, {0.,.5,0.,0.}, {0.,0.,1.,0.} },
            new double[]{1./6,1./3,1./3,1./6},new double[]{1./6,1./3,1./3,1./6});
        if(method.equals("bs23"))return new Tableau(new double[]{0.,.5,.75,1.},
            new double[][]{ {0.,0.,0.,0.}, {.5,0.,0.,0.}, {0.,.75,0.,0.}, {2./9,1./3,4./9,0.} },
            new double[]{2./9,1./3,4./9,0.},new double[]{7./24,1./4,1./3,1./8});
        throw new IllegalArgumentException("unknown method");
    }
    static boolean timesOk(double t,double end,String method) {
        if(!Double.isFinite(end) || end==t)return false;
        double h=end-t,lo=Math.min(t,end),hi=Math.max(t,end),mid=t+.5*h;
        if((method.equals("rk4") || method.equals("bs23")) && !(lo<mid && mid<hi))return false;
        if(method.equals("bs23")) {
            double q=t+.75*h;if(!(Math.min(mid,end)<q && q<Math.max(mid,end)))return false;
        }
        return true;
    }
    static class Trial {
        double[] high,low;int evaluations;
        Trial(double[] high,double[] low,int evaluations){this.high=high;this.low=low;this.evaluations=evaluations;}
    }
    static Trial rkStep(Fn f,double t,double[] y,double end,String method) {
        Tableau tab=tableau(method);double h=end-t;int d=y.length,s=tab.c.length;double[][] stages=new double[s][];
        for(int k=0;k<s;++k) {
            double[] state=new double[d];
            for(int j=0;j<d;++j) {
                double total=0;for(int i=0;i<k;++i)total+=tab.a[k][i]*stages[i][j];state[j]=finite(y[j]+h*total);
            }
            double time=tab.c[k]==1. ? end : finite(t+tab.c[k]*h);double[] slope=f.apply(time,state.clone());
            if(slope.length!=d)throw new IllegalArgumentException("right-hand-side dimension mismatch");
            stages[k]=slope.clone();for(double v:stages[k])finite(v);
        }
        double[] high=new double[d],low=new double[d];
        for(int j=0;j<d;++j) {
            double s1=0,s2=0;for(int k=0;k<s;++k){s1+=tab.b[k]*stages[k][j];s2+=tab.low[k]*stages[k][j];}
            high[j]=finite(y[j]+h*s1);low[j]=finite(y[j]+h*s2);
        }
        return new Trial(high,low,s);
    }
    static double[] pathRow(double t,double[] y) {
        double[] row=new double[y.length+1];row[0]=t;System.arraycopy(y,0,row,1,y.length);return row;
    }
    static class Result {
        double t;double[] y;ArrayList<double[]> path;int accepted,rejected;long evaluations;String status;
        Result(double t,double[] y,ArrayList<double[]> path,int accepted,int rejected,long evaluations,String status) {
            this.t=t;this.y=y;this.path=path;this.accepted=accepted;this.rejected=rejected;this.evaluations=evaluations;this.status=status;
        }
    }
    static Result fixed(Fn f,double t0,double t1,double[] y0,int n,String method) {
        double[] y=validate(t0,t1,y0);
        if((!method.equals("euler") && !method.equals("heun") && !method.equals("rk4")) || n<1)
            throw new IllegalArgumentException("valid fixed method and positive step count required");
        double t=t0;ArrayList<double[]> path=new ArrayList<>();path.add(pathRow(t,y));long calls=0;
        if(t0==t1)return new Result(t,y,path,0,0,0,"done");
        for(int i=0;i<n;++i) {
            double end=i==n-1 ? t1 : finite(t0+(t1-t0)*((double)(i+1)/n));
            if(!timesOk(t,end,method))return new Result(t,y,path,i,0,calls,"stalled");
            Trial trial=rkStep(f,t,y,end,method);y=trial.high;calls+=trial.evaluations;t=end;path.add(pathRow(t,y));
        }
        return new Result(t,y,path,n,0,calls,"done");
    }
    static Result adaptive(Fn f,double t0,double t1,double[] y0,double h0,double atol,double rtol,int limit) {
        double[] y=validate(t0,t1,y0);
        if(!Double.isFinite(h0) || h0<=0 || !Double.isFinite(atol) || atol<=0 || !Double.isFinite(rtol) || rtol<0 || limit<0)
            throw new IllegalArgumentException("invalid initial step, tolerances, or attempt budget");
        double t=t0;ArrayList<double[]> path=new ArrayList<>();path.add(pathRow(t,y));int accepted=0,rejected=0;long calls=0;
        if(t0==t1)return new Result(t,y,path,0,0,0,"done");
        double direction=t1>t0 ? 1. : -1.,size=h0;
        for(int attempt=0;attempt<limit;++attempt) {
            double remaining=Math.abs(t1-t);size=Math.min(size,remaining);
            double end=size==remaining ? t1 : finite(t+direction*size);
            if(!timesOk(t,end,"bs23"))return new Result(t,y,path,accepted,rejected,calls,"stalled");
            Trial trial=rkStep(f,t,y,end,"bs23");calls+=trial.evaluations;double error=0;
            for(int j=0;j<y.length;++j) {
                double scale=finite(atol+rtol*Math.max(Math.abs(y[j]),Math.abs(trial.high[j])));
                error=Math.max(error,Math.abs(finite(trial.high[j]-trial.low[j]))/scale);
            }
            double h=Math.abs(end-t);
            if(error<=1.) {
                t=end;y=trial.high;++accepted;path.add(pathRow(t,y));
                if(t==t1)return new Result(t,y,path,accepted,rejected,calls,"done");
            } else ++rejected;
            double factor=error==0 ? 5. : Math.max(.2,Math.min(5.,.9*Math.pow(error,-1./3)));
            if(error>1.)factor=Math.min(1.,factor);
            size=Math.min(h*factor,Math.abs(t1-t));
        }
        return new Result(t,y,path,accepted,rejected,calls,"limit");
    }
    public static void main(String[] args) {
        Locale.setDefault(Locale.ROOT);Fn f=(t,y) -> new double[]{-y[0]};double exact=Math.exp(-1.);
        ArrayList<Result> results=new ArrayList<>();
        for(String method:new String[]{"euler","heun","rk4"})results.add(fixed(f,0,1,new double[]{1.},4,method));
        results.add(adaptive(f,0,1,new double[]{1.},.5,1e-7,1e-6,10000));
        String[] names={"euler","heun","rk4","adaptive-bs23"};double[] bounds={.06,.005,2e-5,2e-6};boolean ok=true;
        for(int i=0;i<4;++i) {
            Result r=results.get(i);double error=Math.abs(r.y[0]-exact);ok=ok && r.status.equals("done") && r.t==1. && error<=bounds[i];
            System.out.printf("%s: y(1) = %.6f, actual error = %.6f%n",names[i],r.y[0],error);
        }
        Result r=results.get(3);System.out.println("adaptive status = "+r.status);
        System.out.printf("adaptive accepted = %d, rejected = %d, evaluations = %d%n",r.accepted,r.rejected,r.evaluations);
        System.out.println("check = "+(ok?"PASS":"FAIL"));
    }
}

The identical sample output is:

• euler: y(1) = 0.316406, actual error = 0.051473
• heun: y(1) = 0.372529, actual error = 0.004650
• rk4: y(1) = 0.367894, actual error = 0.000015
• adaptive-bs23: y(1) = 0.367879, actual error = 0.000001
• adaptive status = done
• adaptive accepted = 29, rejected = 2, evaluations = 124
• check = PASS

The driver checks method-specific actual-error thresholds and requires every run to reach the target with done. Its coarse Euler result is visibly less accurate; the shared success label means each demonstration met its stated threshold. Function-call count matters as well as accepted steps: the adaptive run used $4(29+2)=124$ evaluations, including rejected work. Other equations or library controllers need not produce those counts.

Complexity, Correctness, and Invariants

Let $d$ be the state dimension, $C_f$ the cost of one vector right-hand-side evaluation, $N$ the fixed step count, $A$ the adaptive attempt count, $K$ its attempt budget, and $S$ the number of accepted steps. Orders below assume sufficient smoothness and a stable convergent regime on a fixed finite interval. The local-error column uses the unnormalized one-step error, not that error divided by $h$.

For a valid nonzero fixed run with representable stages, best, average, and worst loop work at a fixed $N$ are the same. A stalled grid can return earlier. All solvers return immediately without evaluations on a validated zero-length interval. An adaptive nonzero interval can finish in one accepted attempt, but its worst budgeted work includes $K$ unsuccessful or successful attempts; it may return limit. There is no distribution-free average attempt count or unconditional finite count to achieve a tolerance. Stiffness, singularities, rapidly varying derivatives, and evaluator noise can make acceptance difficult.

Consistency and convergence. Euler follows the first Taylor term. Heun and RK4 choose stage locations and weights that match the exact solution's Taylor expansion through higher orders. Under suitable regularity and a Lipschitz stability estimate, errors introduced at each step accumulate to their stated global orders. On well-behaved fixed-interval tests, halving a sufficiently small fixed step reduces global error by factors approaching $2$, $4$, and $16$, respectively. These are asymptotic statements, not universal ratios at coarse steps.

Stage invariant. Every trial starts from a fixed committed $(t,\mathbf y)$. Each explicit stage uses only that original state and already computed stage slopes. The candidate weighted sum does not modify the original state while stages are being evaluated. The tableau arrays in the code encode precisely those dependencies.

Acceptance invariant. The trajectory begins with the initial condition and contains exactly one additional row per accepted step. Its times advance in the requested direction, its last row equals the reported committed state, and a rejected trial never appears in it. Every adaptive accepted trial satisfies its computed normalized-error test. This is a bookkeeping and local diagnostic invariant; it is not a proof of a requested global solution-error bound.

Stability is separate from accuracy order. On $y'=\lambda y$, Euler gives $y_{n+1}=(1+h\lambda)y_n$. With negative real $\lambda$, decay requires $\lvert1+h\lambda\rvert\lt1$. Heun and RK4 have different stability polynomials, but their explicit stability regions are still bounded. Reducing a local error diagnostic does not remove stability restrictions or guarantee that a stiff problem is being solved efficiently.

Variants, Pitfalls, and Optimizations

A stable equation can have an unstable numerical trajectory. For $y'=-10y$ and $y(0)=1$, the exact solution decays. Euler with $h=0.3$ has multiplier $-2$, producing $1,-2,4$ at times $0,0.3,0.6$. Euler with $h=0.05$ has multiplier $1/2$ and decays instead. The smaller step is stable for this test, although stability alone does not make its decay quantitatively exact.

For this equation Euler's strict decay interval is $0\lt h\lt0.2$. At $h=0.2$ the multiplier is $-1$, so the numerical solution oscillates without damping. In genuinely stiff systems, slow modes may coexist with much faster decaying modes, forcing explicit methods to take far smaller steps than the behavior of interest suggests. Implicit Euler, implicit Runge–Kutta methods, Radau, or BDF methods can be appropriate alternatives; they require solving algebraic systems and have their own accuracy and stability contracts.

Avoid the common stage bugs. Do not multiply by $h$ twice when stage variables are derivatives. Do not compute Heun's second slope at the old state, reuse RK4's second stage as its third, or update vector components sequentially in place. For a coupled system such as an oscillator, every component of a stage must use the same stage vector. A scalar example alone may fail to expose these errors.

Treat time and tolerances carefully. Shorten the final trial rather than overshooting the endpoint, and use the actual clipped step in every stage and error calculation. Respect backward-integration direction. Very small time steps can collapse to the current time or make interior stage times coincide; the code reports stalled. A finite accepted state does not imply positivity, energy conservation, or another physical invariant. Scalar absolute tolerances may be unsuitable for state components measured in different units.

Know what adaptive control omits. The maximum norm here differs from the weighted RMS norm used by some libraries. The supplied driver has no dense-output interpolant, event locator, discontinuity handling, or automatic stiffness detection. If a rate law changes at a known time, end a step there and restart with the appropriate model. If a simulation stops on a state-dependent event, bracket and locate that event rather than simply taking the first step past it. Finite-time blow-up can exhaust the step budget or precision; it is not repaired by hiding the failure status.

Embedded estimates can agree accidentally. For $y'=-2y$, $y(0)=1$, a single trial with $h=1/2$ gives both Bogacki–Shampine candidates equal to $1/3$ in exact arithmetic. The exact endpoint is $e^{-1}\approx0.367879$, so the actual error is about $0.034546$ despite a zero embedded difference. On the linear test equation, their difference is $-y z^3(1+z)/48$ with $z=h\lambda$, which vanishes at $z=-1$. A tight tolerance cannot detect that particular coincidence; compare smaller initial steps or an independent method when validating a solve. This is another reason the acceptance test is a diagnostic rather than a global certificate.

Choose a method for the problem. Dormand–Prince $5(4)$ is another common adaptive embedded pair; its stage coefficients and controller order differ from the $3(2)$ pair implemented here. Symplectic integrators can better preserve long-term geometric behavior in Hamiltonian systems, while general-purpose RK4 does not conserve energy exactly. Higher-order equations should be converted to first-order systems before using these interfaces.

Improve work without changing the method. Bogacki–Shampine's final stage is the derivative at the accepted endpoint, so a careful implementation can reuse it as the next first stage. Cache or batch other expensive shared calculations when the rate law permits. If only the final state is needed, omit trajectory storage to reduce memory to $O(d)$. Keep accepted and rejected evaluation costs visible when comparing solvers.

Practice Problems and Real-World Use

Start with exact solutions and step-refinement experiments, then add coupled states, rejection tests, and difficult dynamics. Record actual error when known, accepted and rejected counts, evaluation counts, endpoint time, and any invariant relevant to the model.

For further reading, the authors' Fundamentals of Numerical Computation chapter on Runge–Kutta methods explains stage construction and order. SciPy's RK23 documentation identifies the Bogacki–Shampine $3(2)$ pair, and its initial-value solver documentation distinguishes explicit and implicit solver choices. Library norms, derivative reuse, interpolation, and stopping controls can differ from this explicit teaching implementation.