Numerical Integration of Differential Equations in MATLAB: Benchmarking Accuracy, Stability, Stiffness, and Conservation

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, conservation laws, and long-time orbital dynamics. This article presents a comprehensive discussion of differential equation integrator benchmarking in MATLAB. The discussion is based on the project: GitHub Repository: https://github.com/mohammadijoo/DifferentialEquationIntegratorBenchmarks_MATLAB The repository implements and compares a broad collection of numerical ODE integrators on classical benchmark problems. A companion YouTube simulation video demonstrates the MATLAB implementation, generated plots, MATLAB source codes, and exported CSV metric files. Video: MATLAB implementation of differential equation integrator benchmarks, including repository code execution, generated plots, and CSV metric files. 1. Initial-Value Problems and the Role of Numerical Integrators The main mathematical object in this benchmark framework is the initial-value problem: d t d y ​ = f ( t , y ) , y ( t 0 ​ ) = y 0 ​ Here, the state vector belongs to an n-dimensional real space: y ( t ) ∈ R n The function 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: t n ​ = t 0 ​ + nh where h is the time step. The numerical approximation is denoted by: y n ​ ≈ y ( t n ​ ) The purpose of a benchmark is to compare different update rules of the form: y n + 1 ​ = Φ h ​ ( t n ​ , y n ​ ) where the map Phi_h is the numerical flow map produced by a specific method. A good 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 remain efficient? Does it fail, blow up, or produce NaN values? Does it preserve qualitative dynamics such as oscillations, attractors, or orbits? These questions are especially important in MATLAB, where numerical simulation is frequently used in control engineering, robotics, mechanical dynamics, electrical circuits, chemical kinetics, and computational physics. 2. Benchmarking Is Not Just About Error The most direct error metric is the final-time error: e final ​ = ∣ y N ​ − y ref ​ ( T ) ∣ where y_ref(T) is either the exact solution or a high-accuracy reference solution at the final time T . However, this metric alone can be misleading. A method may have small final error but large intermediate error. Therefore, a trajectory-level maximum error is also useful: e m a x ​ = 0 ≤ n ≤ N max ​ ∣ y n ​ − y ref ​ ( t n ​ ) ∣ A root-mean-square error gives a global average measure: e RMSE ​ N + 1 1 ​ n = 0 ∑ N ​ ∣ y n ​ − y ref ​ ( t n ​ ) ∣ 2 ​ For computational efficiency, the benchmark should also record CPU time: T CPU ​ number of right-hand-side evaluations: N fev ​ number of accepted steps: N steps ​ number of rejected adaptive steps: N rejected ​ number of Newton iterations for implicit methods: N Newton ​ and number of Jacobian evaluations: N Jacobian ​ For conservative systems, an additional invariant error is critical. If a