Numerical integration is one of the central tools of scientific computing. Whenever a physical, biological, chemical, mechanical, electrical, or control system is modeled by differential equations, the next practical question is usually not only:

What is the equation?

but also:

Which numerical method should be trusted to simulate it?

This question is more subtle than it first appears.

A numerical integrator that performs very well on a smooth non-stiff control example may fail completely on a stiff chemical-kinetics model. A method that gives excellent short-time accuracy may slowly destroy the energy of a mechanical system. A method that is computationally cheap may require such a small time step that it becomes inefficient. A method that is stable for one equation may be unstable for another.

For this reason, numerical integrators should be studied through benchmark equations that expose different mathematical behaviors:

• linear stability,
• oscillation,
• nonlinear stiffness,
• chaos,
• chemical reaction stiffness,
• conservation laws,
• long-time Hamiltonian behavior,
• positivity constraints,
• PDE stiffness after spatial discretization,
• and failure modes such as blow-up or non-applicability.

This article is based on my MATLAB repository:

GitHub Repository:
https://github.com/mohammadijoo/DifferentialEquationIntegratorBenchmarks_MATLAB

The repository started as a MATLAB ODE integrator benchmark suite, but it has now been reorganized into a broader numerical-solver catalog. It includes runnable MATLAB benchmark code, implemented integrators, benchmark documentation, generated results, and a much larger roadmap for future DAE, PDE, SDE, DDE, fractional, BVP, hybrid, multirate, and parallel-in-time benchmark directions.

1. The Main Mathematical Object: Initial-Value Problems

Most of the currently runnable benchmark framework is based on initial-value problems of the form:

where

is the state vector, f(t,y) is the right-hand-side function, and y_0 is the initial condition.

A numerical integrator approximates the exact solution at discrete time points:

and produces numerical states

A one-step numerical method can be written abstractly as

where Phi_h is the numerical flow map generated by the method.

A benchmark should not only ask whether y_n is close to the exact solution. It should also ask:

• Is the method stable?
• Does it preserve physical invariants?
• Does it handle stiffness?
• Does it preserve positivity?
• Does it require many function evaluations?
• Does it require many Newton or linear-solve iterations?
• Does it fail, blow up, or return non-applicable?
• Does it preserve qualitative behavior such as oscillations, orbits, phase portraits, or attractors?

These questions are especially important in MATLAB, where numerical simulation is widely used in control engineering, robotics, mechanical dynamics, electrical circuits, chemical kinetics, computational physics, and scientific computing education.

2. What Changed in the Repository

The original version of the project focused mainly on a smaller set of classical ODE benchmarks and common integrators.

The current repository is broader. It now has two related layers.

Layer 1: Runnable MATLAB Benchmark Framework

This is the part that can be executed directly from MATLAB.

The main entry points are:

run_all

and

run_one_benchmark('LinearTestEquation')
run_one_benchmark('LorenzSystem')
run_one_benchmark('RobertsonKinetics')
run_one_benchmark('KeplerTwoBody')

The default options are obtained from:

opts = default_options();

For example:

opts = default_options();
opts.h = 1e-2;
opts.write_plots = true;
opts.write_tables = true;
opts.include_planned_methods = false;
opts.include_planned_benchmarks = false;

results = run_all(opts);

By default, planned methods and planned benchmarks are excluded:

opts.include_planned_methods = false;
opts.include_planned_benchmarks = false;

This is important. The repository contains many documented and scaffolded future methods, but the default runner only executes methods and benchmarks marked as implemented.

Layer 2: Expanded Documentation and Roadmap

The repository has also become a structured catalog of numerical methods and benchmark problems.

The documentation and source folders are organized by method and benchmark families. The method catalog includes:

• classical ODE initial-value integrators,
• stiff, implicit, multistep, Rosenbrock, exponential, IMEX, splitting, and geometric methods,
• DAE integrators and constraint methods,
• PDE time integration and spatial discretization methods,
• stochastic differential equation methods,
• delay differential equation methods,
• fractional differential equation methods,
• boundary-value problem methods,
• multirate and parallel-in-time methods.

This means the repository is not only a runnable MATLAB benchmark suite, but also an educational numerical-analysis roadmap.

3. Benchmarking Is Not Just About Final-Time Error

The most direct metric is the final-time error:

where y_ref(T) is either an exact solution or a high-accuracy reference solution at the final time T.

However, final-time error alone can be misleading. A method may accidentally have small error at the final time but large error during the trajectory. Therefore, a trajectory-level maximum error is also useful:

A root-mean-square trajectory error gives an average global measure:

For computational cost, the benchmark records or emphasizes quantities such as:

• CPU time,
• number of accepted steps,
• number of right-hand-side evaluations,
• number of rejected adaptive steps,
• number of Jacobian evaluations,
• number of Newton iterations,
• number of linear solves,
• method status.

For implicit and stiff methods, CPU time alone is not enough. A solver with more steps may still be fast if each step is cheap. A solver with fewer steps may be expensive if each step requires Jacobians, Newton iterations, and matrix factorizations.

For conservative or mechanical systems, invariant drift is also essential. If a system has an invariant I(y), such as energy, mass, or angular momentum, then the invariant error can be measured as:

and the maximum invariant drift is

This metric is especially important for oscillators, Hamiltonian systems, celestial mechanics, robotics simulation, and long-time mechanical dynamics.

4. Stability and the Dahlquist Test Equation

The linear scalar test equation is one of the simplest ways to study basic stability behavior:

The exact solution is

If lambda < 0, the exact solution decays. A numerical method should reproduce that decay without artificial growth.

For a one-step method applied to the linear test equation, the numerical update can often be written as

where R(z) is the stability function of the method.

The method is stable for a given z if

This benchmark is useful because it exposes the absolute-stability region of a method. For example, Forward Euler is simple, but its stability region is limited. Backward Euler is much more robust for decay problems and is often useful for stiff systems, but it can strongly damp dynamics.

5. The Original Six Core Benchmarks

The early version of the article focused on six major ODE benchmarks. These are still important and remain part of the project.

5.1 Linear Test Equation

Purpose:

• stability behavior,
• decay accuracy,
• basic convergence,
• step-size sensitivity.

Equation:

This benchmark is small, but it is not trivial. It is the standard first test for understanding whether a method behaves correctly on exponential growth or decay.

5.2 Harmonic Oscillator

The harmonic oscillator is a basic conservative mechanical system:

Its energy is

For the exact solution,

This benchmark tests:

• phase accuracy,
• amplitude error,
• energy drift,
• long-time oscillatory behavior,
• suitability of symplectic or energy-preserving methods.

A method can have small local error but still slowly increase or decrease the oscillator energy. That is why geometric diagnostics are important.

5.3 Van der Pol Oscillator

The Van der Pol oscillator is a nonlinear oscillator:

For small mu, the problem behaves like a moderately nonlinear oscillator. For large mu, it develops fast and slow time scales and becomes stiff.

This benchmark tests:

• nonlinear dynamics,
• limit-cycle behavior,
• stiffness transition,
• robustness under changing time scales.

It is useful because a method that works well for mu = 1 may struggle when mu becomes large.

5.4 Lorenz System

The Lorenz system is a classical chaotic system:

This benchmark tests:

• nonlinear instability,
• sensitive dependence on initial conditions,
• qualitative trajectory behavior,
• attractor visualization,
• reliability over finite time windows.

For chaotic systems, a numerical trajectory should not be expected to match a reference solution forever. Instead, the benchmark should focus on reasonable finite-time accuracy and qualitative behavior.

5.5 Robertson Stiff Chemical Kinetics

The Robertson problem is a classical stiff chemical kinetics benchmark:

It also has a mass-conservation property:

and the species concentrations should remain nonnegative:

This benchmark tests:

• stiffness,
• positivity preservation,
• conservation of total mass,
• robustness of implicit and linearly implicit methods,
• cost of Jacobians, Newton iterations, and linear solves.

Explicit methods can require extremely small time steps on this problem. Stiff solvers are usually more appropriate.

5.6 Kepler Two-Body Problem

The Kepler two-body problem describes motion under an inverse-square central force:

where

The total energy is

The angular momentum in two dimensions is

For the exact Kepler problem,

This benchmark tests:

• long-time orbital behavior,
• energy drift,
• angular momentum drift,
• phase error,
• orbit distortion,
• suitability of symplectic methods.

A method may have good short-time accuracy but slowly distort the orbit. In orbital mechanics and robotics simulation, long-time qualitative fidelity is often more important than only pointwise short-time error.

6. Expanded Implemented Benchmark Families

The current repository goes beyond the original six examples. It now includes a wider set of implemented benchmarks grouped by mathematical behavior.

6.1 Non-Stiff Accuracy and Nonlinear Dynamics

This group includes smooth and nonlinear ODE examples such as:

• Brusselator ODE,
• Duffing oscillator,
• forced damped pendulum,
• logistic growth,
• Lotka-Volterra predator-prey model,
• rigid-body Euler equations,
• SIR epidemic model.

These benchmarks are useful for testing:

• nonlinear behavior,
• oscillations,
• phase-space geometry,
• population dynamics,
• sensitivity to step size,
• qualitative trajectory preservation.

For example, the logistic equation has the form:

This benchmark is useful because it has simple nonlinear saturation and a known qualitative behavior: solutions should approach the carrying capacity K.

The Lotka-Volterra predator-prey model has the form:

It tests whether a method can represent coupled nonlinear oscillations without artificial damping or growth.

6.2 Stiff Multiscale and Chemical Kinetics Benchmarks

The stiff benchmark family now includes examples such as:

• Akzo Nobel chemical kinetics,
• Allen-Cahn method-of-lines stiff ODE,
• HIRES problem,
• Oregonator / Belousov-Zhabotinsky model,
• pollution model ODE,
• Prothero-Robinson problem,
• ring modulator circuit ODE,
• stiff Brusselator regime.

These examples are important because stiffness is not only a theoretical issue. It appears in chemical kinetics, reaction-diffusion systems, circuit models, and multiscale physical systems.

A stiff system usually contains widely separated time scales. Informally, a system may have both fast decaying modes and slow solution features. The numerical method must remain stable on the fast modes while accurately tracking the slow dynamics.

A linearized stiff system may have eigenvalues satisfying

where the fastest stable mode imposes a severe step-size restriction on explicit methods.

This is why stiff benchmarks should record more than error and CPU time. They should also examine:

• rejected steps,
• Jacobian evaluations,
• Newton iterations,
• linear solves,
• failure status,
• positivity violation,
• conservation error.

6.3 Hamiltonian, Geometric, and Long-Time Mechanics Benchmarks

The repository also includes mechanical and Hamiltonian-style benchmarks such as:

• charged particle in magnetic field,
• double pendulum,
• FPUT lattice,
• Henon-Heiles Hamiltonian,
• N-body gravitational problem,
• planar restricted three-body problem,
• Toda lattice.

These benchmarks are designed to expose long-time structure preservation.

A Hamiltonian system has the canonical form:

For the exact continuous system, the Hamiltonian is conserved:

However, a general-purpose integrator does not necessarily preserve this structure. A non-symplectic method may introduce artificial damping or artificial energy growth over long simulations.

This is why symplectic methods, variational integrators, discrete-gradient methods, and energy-momentum ideas are important for this family of benchmarks.

6.4 PDE Method-of-Lines and Conservation-Law Benchmarks

The expanded repository also includes PDE-derived problems after spatial discretization, including:

• advection equation,
• Cahn-Hilliard equation,
• Fisher-KPP equation,
• Gray-Scott reaction-diffusion,
• heat equation,
• inviscid Burgers equation,
• Korteweg-de Vries equation,
• Kuramoto-Sivashinsky equation,
• viscous Burgers equation,
• wave equation.

A PDE can often be converted into a large ODE system using the method of lines. A general semi-discrete form is:

or, when the mass matrix is the identity,

For diffusion-dominated problems, spatial discretization can generate stiffness. For hyperbolic conservation laws, numerical dissipation, dispersion, shock handling, and stability restrictions become important.

For example, the heat equation

after finite-difference spatial discretization becomes a system of ODEs. The eigenvalues of the discrete Laplacian can impose severe time-step restrictions on explicit methods.

The advection equation

tests phase accuracy, numerical diffusion, and wave propagation behavior.

These PDE-derived benchmarks make the project more useful for computational physics, numerical PDEs, and scientific machine simulation workflows.

7. Planned Benchmark Directions

The repository also contains planned or scaffolded benchmark directions. These should not be treated as fully validated runnable benchmark cases until their implementation status is changed to implemented.

The planned directions include:

• DAE constrained systems and circuits,
• stochastic differential equations,
• delay differential equations,
• fractional differential equations,
• boundary-value problems,
• hybrid and event-driven systems.

This planned layer is useful because it shows how the project can grow beyond ordinary ODE integration.

For a DAE, the problem may have the form:

For a delay differential equation:

For an SDE:

For a fractional differential equation with a Caputo derivative:

Each of these problem classes requires different diagnostics. For example, SDEs require statistical metrics, DDEs require history interpolation diagnostics, DAEs require constraint residuals, and fractional equations require memory-aware numerical methods.

8. Implemented Integrator Families

The current runnable method registry includes several important families of numerical methods.

8.1 Explicit One-Step and Low-Order Runge-Kutta Methods

Implemented examples include:

• Forward Euler,
• Explicit midpoint,
• Heun method,
• Ralston method,
• Kutta RK3,
• classical RK4.

Forward Euler is the simplest method:

It is first order:

It is useful as a baseline, but it is usually not competitive for serious simulation because of its low order and limited stability.

The explicit midpoint method uses an internal half-step:

Classical RK4 uses four stages:

RK4 is a strong baseline for smooth non-stiff ODEs, but it is not a stiff solver and not a structure-preserving geometric method.

8.2 Adaptive Runge-Kutta Methods

Implemented adaptive RK examples include:

• Runge-Kutta-Fehlberg RKF45,
• Dormand-Prince RK45.

Adaptive methods estimate local error by comparing two approximations of different orders. A typical embedded pair produces

and estimates the local error as

The step size is then adjusted. A common form is

where s is a safety factor.

Adaptive RK methods are useful for non-stiff problems with varying solution scales. However, adaptivity does not automatically solve stiffness. A stiff problem may still require implicit, BDF, Rosenbrock, or other stiffly stable methods.

8.3 Implicit One-Step Methods

Implemented implicit one-step examples include:

• Backward Euler,
• implicit midpoint,
• implicit trapezoidal rule.

Backward Euler is

This is an implicit equation for y_{n+1}. It usually requires a nonlinear solve.

Define the residual:

The method solves

typically by Newton's method.

Backward Euler is first order and strongly damping, but it is very useful for stiff decay problems.

The implicit trapezoidal rule is

It is second order and often less dissipative than Backward Euler.

8.4 Linear Multistep and Predictor-Corrector Methods

Implemented multistep examples include:

• Adams-Bashforth 2,
• Adams-Bashforth 4,
• Adams-Moulton 2,
• Adams-Bashforth-Moulton 4 predictor-corrector,
• BDF2.

A linear multistep method uses multiple previous solution values:

Adams-Bashforth methods are explicit. BDF methods are implicit and are important in stiff ODE integration.

BDF2 has the form

Because BDF2 is implicit, it generally requires nonlinear solves, Jacobians, and linear algebra. Therefore, benchmarking BDF-style methods should include more than error and CPU time.

8.5 Rosenbrock and Linearly Implicit Methods

Implemented linearly implicit examples include:

• Rosenbrock-Euler,
• ROS2.

Rosenbrock methods avoid fully nonlinear solves by linearizing around the current state. A simplified Rosenbrock-Euler form is:

where

These methods are useful for stiff systems because they use Jacobian information without requiring a full nonlinear Newton iteration at every step.

The benchmark should record:

• Jacobian evaluations,
• linear solves,
• right-hand-side evaluations,
• stability behavior,
• and failure status.

8.6 Exponential Methods

The implemented catalog also includes exponential Euler.

Exponential methods are important for systems with a stiff linear part:

The exact solution formula involves the matrix exponential:

Exponential Euler approximates the nonlinear part while treating the linear part through exponential propagation.

These methods are especially relevant for semi-discretized PDEs, where the linear stiff part may dominate the stability restriction.

8.7 Geometric, Symplectic, and Energy-Aware Methods

Implemented geometric methods include:

• Velocity Verlet,
• Stormer-Verlet,
• Yoshida fourth-order symplectic composition,
• discrete-gradient oscillator method.

For a separable Hamiltonian system

the equations are

Velocity Verlet has the form

Symplectic methods are not designed to minimize final-time error only. Their strength is long-time qualitative behavior. They often keep the numerical trajectory on a more physically meaningful path over long simulations.

9. Work-Precision Behavior

A useful benchmark should compare error against cost.

A work-precision diagram usually plots an error metric against a work metric such as:

• CPU time,
• number of function evaluations,
• number of accepted steps,
• number of Jacobian evaluations,
• number of linear solves.

A method is not better simply because it has lower error at one step size. It should be compared across different tolerances or step sizes.

A typical work-precision study asks:

• How fast does error decrease as work increases?
• Which methods fail before reaching the target accuracy?
• Which methods become unstable for larger step sizes?
• Which methods are accurate but too expensive?
• Which methods are cheap but unreliable?

For example, RK4 may be efficient for smooth non-stiff ODEs, but a stiff chemical kinetics problem may require a stiff method even if each step is more expensive.

10. Applicability Matters

One important improvement in the repository is the explicit idea of method applicability.

Not every method should run on every benchmark.

For example:

• Velocity Verlet is meaningful for separable mechanical systems, but not for arbitrary chemical kinetics.
• Stochastic methods are meaningful for SDEs, not deterministic ODEs.
• DAE methods are meaningful for constrained problems, not ordinary explicit ODE-only tests.
• PDE spatial discretization methods belong to PDE-derived benchmarks.
• BVP methods should not be evaluated like initial-value solvers.

Therefore, the benchmark framework should distinguish between:

• success,
• failure,
• skipped,
• not applicable,
• planned but not implemented.

This is better than forcing every method to run on every problem. A solver comparison should respect mathematical structure.

11. Repository Structure

The repository now mirrors documentation and source-code organization.

A simplified view is:

.
├── README.md
├── LICENSE
├── CONTRIBUTING.md
├── CITATION.cff
├── run_all.m
├── run_one_benchmark.m
├── docs/
│ ├── methods/
│ ├── benchmarks/
│ ├── tutorials/
│ └── references.md
├── src/
│ ├── config/
│ │ ├── default_options.m
│ │ ├── method_registry.m
│ │ └── benchmark_registry.m
│ ├── core/
│ │ ├── run_benchmark.m
│ │ ├── compute_metrics.m
│ │ └── newton_solve.m
│ ├── methods/
│ ├── benchmarks/
│ └── plotting/
├── tests/
└── results/

The registries are important because they define which methods and benchmarks are implemented, planned, applicable, or skipped.

To inspect implemented methods and benchmarks inside MATLAB:

methods = method_registry();
benchmarks = benchmark_registry();

implementedMethods = methods([methods.implemented]);
implementedBenchmarks = benchmarks([benchmarks.implemented]);

{implementedMethods.name}'
{implementedBenchmarks.name}'

12. Recommended MATLAB Workflow

A typical workflow is:

run_all

To run a single benchmark:

run_one_benchmark('VanDerPolOscillator')

To customize options:

opts = default_options();

opts.h = 5e-3;
opts.n_output = 2001;
opts.write_plots = true;
opts.write_tables = true;
opts.write_global_summary = true;

results = run_one_benchmark('LorenzSystem', opts);

To keep the run restricted to implemented entries:

opts.include_planned_methods = false;
opts.include_planned_benchmarks = false;

Results are written under:

results/

or under benchmark-specific folders such as:

results/<BenchmarkName>/

Typical outputs include:

• CSV summary tables,
• final-error plots,
• CPU-time plots,
• function-evaluation plots,
• work-precision diagrams,
• invariant drift plots,
• phase portraits,
• orbit plots,
• attractor projections,
• benchmark-specific solution visualizations.

13. Why the Expanded Benchmark Set Is Useful

The expanded repository is useful because numerical integrators have different strengths.

Explicit Runge-Kutta methods are often strong for smooth non-stiff ODEs.

Adaptive RK methods are useful when local-error control is needed.

Implicit methods are important for stiff decay and stiff chemical kinetics.

BDF and multistep methods are important for stiff and long-time integration.

Rosenbrock methods are useful when Jacobian-based linearized stiffness handling is needed.

Exponential methods are useful when the problem has a strong linear stiff part.

Symplectic and geometric methods are important for long-time Hamiltonian mechanics.

PDE method-of-lines benchmarks expose stiffness, dispersion, diffusion, and conservation issues caused by spatial discretization.

DAE, SDE, DDE, fractional, BVP, and hybrid systems require additional diagnostics beyond ordinary ODE error metrics.

Therefore, the main lesson is:

There is no universal best integrator. There is only the best integrator for a given problem structure, accuracy requirement, stability constraint, and computational budget.

14. Practical Interpretation of Benchmark Results

When interpreting the benchmark outputs, I would avoid ranking methods using a single metric.

A better interpretation is problem-dependent.

For the linear test equation, stability behavior and decay accuracy matter.

For the harmonic oscillator, phase error and energy drift matter.

For Van der Pol, stiffness transition and nonlinear oscillation matter.

For Lorenz, finite-time accuracy and qualitative chaotic behavior matter.

For Robertson kinetics, stiffness, positivity, mass conservation, and solver robustness matter.

For Kepler and other Hamiltonian systems, long-time orbital geometry, energy drift, and angular momentum drift matter.

For PDE method-of-lines benchmarks, both spatial and temporal discretization behavior matter.

For future DAE benchmarks, constraint residuals and consistent initialization will matter.

For future SDE benchmarks, strong and weak convergence should be measured statistically, not only through a single deterministic trajectory.

This is why benchmark results should be interpreted as a numerical profile, not as a single winner list.

15. Conclusion

Numerical integration is not only a matter of implementing formulas. It is a matter of matching the numerical method to the mathematical structure of the differential equation.

The same method can be excellent on one benchmark and poor on another. RK4 may be an excellent general-purpose baseline for smooth non-stiff systems, but it is not a stiff solver and not a geometric energy-preserving method. Backward Euler may handle stiff decay robustly, but it can overdamp conservative systems. Symplectic methods may preserve long-time mechanical structure, but they are not designed for arbitrary stiff chemical kinetics. Adaptive RK methods are powerful, but adaptivity does not automatically solve stiffness.

The updated MATLAB repository reflects this idea by separating runnable implementations from a larger numerical-method roadmap. It includes implemented methods, implemented benchmarks, planned/scaffolded future directions, plotting utilities, result tables, and documentation categories for a broad set of differential-equation problems.

The most important conclusion is:

A good benchmark suite should measure accuracy, cost, stability, stiffness behavior, invariant drift, method applicability, and failure modes. Numerical integration is problem-dependent, and solver choice should be guided by the mathematical structure of the model.