# Time stepping
Time stepping is the numerical technique that advances the solution of a time-dependent differential equation from one time level to the next, producing a sequence of approximate values \( y_n \approx y(t_n) \) at discrete times separated by a step \( h \). It is the standard way to solve initial value problems for [ordinary differential equations](https://www.edgepedia.org/ordinary-differential-equation) (ODEs) and, after spatial discretization, the semi-discrete systems arising from partial differential equations (PDEs).[1](https://unige.ch/~hairer/preprints/pcam-ode.pdf) The output is not a closed-form formula but a table of approximate solution values, each computed from earlier ones by an algebraic update rule.[1](https://unige.ch/~hairer/preprints/pcam-ode.pdf)
| Key fact | Statement |
|---|---|
| What it produces | A sequence \( y_n \) approximating the solution at successive times, computed by recursive updates[1](https://unige.ch/~hairer/preprints/pcam-ode.pdf) |
| Forward Euler | \( y_{n+1} = y_n + h f(t_n, y_n) \); first order, explicit[1](https://unige.ch/~hairer/preprints/pcam-ode.pdf) |
| Classical RK4 | Four stages per step, combined as \( u_{n+1} = u_n + \tfrac{1}{6}(k_1 + 2k_2 + 2k_3 + k_4) \); fourth order[2](https://www.emse.fr/~bonnefoy/Public/MetNum-EMSE.pdf) |
| Implicit workhorses | Backward Euler and Crank–Nicolson (trapezoidal, \( \theta = 1/2 \)) require solving equations at each step but are A-stable[3](https://jschoeberl.github.io/IntroSC/ODEs/firstmethods.html) • [4](https://dumux.org/docs/doxygen/master/benchmark-timestepping-methods.html) |
| BDF methods | Implicit multistep schemes of order \( k \), unstable for \( k > 6 \); the standard choice for stiff problems[1](https://unige.ch/~hairer/preprints/pcam-ode.pdf) |
| Second Dahlquist barrier | An A-stable linear multistep method cannot have order above 2[5](http://www.scholarpedia.org/article/Linear_multistep_method) |
| Adaptive control | An embedded pair estimates the local error and rescales the step by \( \tilde{h} = h \, (\mathrm{tol}/\|\ell_n\|)^{1/(p+1)} \) with a safety factor below 1[6](https://drreynolds.github.io/files/viasm-2025-tutorial.pdf) |
## How it works
Discretizing time turns a differential equation \( y' = f(t, y) \) into a recursion: each step replaces the derivative by algebraic combinations of solution and slope values. Deriving the update from the integral form \( y(t_{n+1}) = y(t_n) + \int_{t_n}^{t_{n+1}} f \, dt \) makes the family structure transparent: the left rectangular rule gives forward Euler, the right rectangular rule gives backward Euler, and the [trapezoidal rule](https://www.edgepedia.org/trapezoidal-rule) gives Crank–Nicolson.[3](https://jschoeberl.github.io/IntroSC/ODEs/firstmethods.html) A method is explicit if the new value appears only through known quantities, and implicit if it appears inside \( f \), so that each step requires solving a (generally nonlinear) system, typically with [Newton's method](https://www.edgepedia.org/newtons-method).[3](https://jschoeberl.github.io/IntroSC/ODEs/firstmethods.html)
Accuracy is measured by order. A method has order \( p \) if the local error of one step from the exact solution is bounded by \( C \cdot h^{p+1} \); the errors of \( O(1/h) \) steps then accumulate to a global error \( O(h^p) \).[1](https://unige.ch/~hairer/preprints/pcam-ode.pdf) Convergence theory closes the loop: Dahlquist's equivalence theorem states that a consistent linear multistep method converges if and only if it is stable, and the Lax–Richtmyer equivalence gives the analogous statement for well-posed linear IVPs.[7](https://bpb-us-w2.wpmucdn.com/sites.brown.edu/dist/1/376/files/2022/04/Chapter-3-Review-of-Time-Stepping-Algorithms-for-ODEs.pdf) • [8](https://djps.github.io/docs/numericalmethods/part3/odes/)
Stability is analyzed on the test equation \( y' = \lambda y \) with \( \mathrm{Re}(\lambda) < 0 \): a method with stability function \( R(z) \) is A-stable if \( |R(z)| \le 1 \) for the entire left half-plane. Forward Euler has \( R(z) = 1 + z \) and is not A-stable; backward Euler has \( R(z) = (1-z)^{-1} \) and is.[9](https://www.csc.univie.ac.at/files/Numerical_methods_for_differential_equations_W2015/main_1.pdf) A scheme is conditionally stable when stability depends on \( h \) (and possibly the spatial mesh), unconditionally stable when any \( h \) is safe.[2](https://www.emse.fr/~bonnefoy/Public/MetNum-EMSE.pdf)
Stiffness is the practical mismatch: the step size needed to keep an explicit method stable is much smaller than the one needed to represent the solution accurately.[6](https://drreynolds.github.io/files/viasm-2025-tutorial.pdf) For the heat equation, explicit Euler forces \( h = O(\Delta x^2) \), a restriction usually considered very restrictive; semi-implicit schemes remove it at the price of solving a large sparse linear system each step.[10](https://link.springer.com/article/10.1186/s40323-016-0063-y) For hyperbolic problems, explicit step sizes are limited by the finest mesh patch or fastest wave speed, the CFL-type constraint.[11](https://excalibur-neptune.github.io/Documents/_static/TN-03_AReviewTimeSteppingTechniquesPreconditioningHyperbolicAnisotropicEllipticProblem.pdf) Barriers limit implicit multistep methods: A-stability forces order \( r \le 2 \),[5](http://www.scholarpedia.org/article/Linear_multistep_method) and [BDF](https://www.edgepedia.org/backward-differentiation-formula) methods of order \( k \le 6 \) are A(α)-stable but become unstable for \( k > 6 \).[5](http://www.scholarpedia.org/article/Linear_multistep_method) No such order barrier exists for A-stable Runge–Kutta methods.[1](https://unige.ch/~hairer/preprints/pcam-ode.pdf)
## How it is done
The θ-family \( y_{n+1} = y_n + h\big[(1-\theta) f(t_n, y_n) + \theta f(t_{n+1}, y_{n+1})\big] \) spans the basic one-step schemes: \( \theta = 0 \) gives explicit Euler, \( \theta = 1 \) backward Euler, and \( \theta = 1/2 \) Crank–Nicolson.[7](https://bpb-us-w2.wpmucdn.com/sites.brown.edu/dist/1/376/files/2022/04/Chapter-3-Review-of-Time-Stepping-Algorithms-for-ODEs.pdf) An \( s \)-stage Runge–Kutta method evaluates \( s \) slopes \( k_i = f(t_n + c_i \cdot h,\; u_n + h \sum_j a_{i,j} \cdot k_j) \); it is explicit exactly when \( a_{i,j} \) is strictly lower triangular, so that \( a_{i,j} = 0 \) for \( j \ge i \).[8](https://djps.github.io/docs/numericalmethods/part3/odes/) The classical fourth-order method uses
\[ k_1 = h f(t_n, u_n), \quad k_2 = h f(t_n + \tfrac{h}{2}, u_n + \tfrac{k_1}{2}), \quad k_3 = h f(t_n + \tfrac{h}{2}, u_n + \tfrac{k_2}{2}), \quad k_4 = h f(t_n + h, u_n + k_3), \]
\[ u_{n+1} = u_n + \tfrac{1}{6}(k_1 + 2 k_2 + 2 k_3 + k_4). \]
It costs four slope evaluations per step.[2](https://www.emse.fr/~bonnefoy/Public/MetNum-EMSE.pdf) Among implicit multistep schemes, BDF2 reads \( (3u^{n+1} - 4u^n + u^{n-1})/(2h) = f(t^{n+1}, u^{n+1}) \), combining second-order accuracy with strong stability.[2](https://www.emse.fr/~bonnefoy/Public/MetNum-EMSE.pdf) Order limits constrain design: an \( s \)-stage explicit RK method cannot exceed order \( s \), and implicit RK order cannot exceed \( 2s \).[8](https://djps.github.io/docs/numericalmethods/part3/odes/)
The choice of scheme follows the stiffness diagnosis. Stiff equations require implicit methods, and stiffly stable BDF schemes are recommended for efficiency; for nonstiff problems, explicit [Runge–Kutta methods](https://www.edgepedia.org/runge-kutta-methods) of orders up to 8 and Adams-type multistep methods up to order 12 are the most widely used.[7](https://bpb-us-w2.wpmucdn.com/sites.brown.edu/dist/1/376/files/2022/04/Chapter-3-Review-of-Time-Stepping-Algorithms-for-ODEs.pdf) • [1](https://unige.ch/~hairer/preprints/pcam-ode.pdf) Because explicit methods are stability limited, practitioners often pick a step size just small enough for linear stability; implicit methods are typically not stability limited.[6](https://drreynolds.github.io/files/viasm-2025-tutorial.pdf)
[Adaptive control](https://www.edgepedia.org/adaptive-control) works by estimating the local truncation error each step, accepting the step if it is small enough and otherwise recomputing with a new step size.[12](https://drreynolds.github.io/files/atpesc-2022-tutorial.pdf) Embedded Runge–Kutta pairs add a second set of \( b \) coefficients that reuse the stored stage vectors \( k_i \) to produce a lower-order solution; the difference \( \|y_{n+1} - \tilde{y}_{n+1}\| \) estimates the error, and the controller updates \( \tilde{h} = h_n (\mathrm{tol}/\|\ell_n\|)^{1/(p+1)} \) with a safety factor below 1.[6](https://drreynolds.github.io/files/viasm-2025-tutorial.pdf) Embedded pairs of orders 5 and 8 in the Dormand–Prince family are standard for this purpose.[1](https://unige.ch/~hairer/preprints/pcam-ode.pdf) For PDEs, the [method of lines](https://www.edgepedia.org/method-of-lines) reduces the problem to an ODE or DAE system, to which general-purpose packages apply: SUNDIALS provides CVODE for IVPs, ARKODE for linearly implicit, split, and multirate (MRIStep) integration, and IDA for DAEs; PETSc's TS module offers a unified implicit, explicit, and IMEX interface.[12](https://drreynolds.github.io/files/atpesc-2022-tutorial.pdf)
## Origin
The development of time stepping rests on a small set of papers. C. F. Curtiss and J. O. Hirschfelder's 1952 paper in the [Proceedings of the National Academy of Sciences](https://www.edgepedia.org/proceedings-of-the-national-academy-of-sciences), "Integration of Stiff Equations," is associated with stiff differential equations and the BDF methods used to solve them.[13](https://doi.org/10.1073/pnas.38.3.235) Germund G. Dahlquist's 1963 BIT Numerical Mathematics paper, "A special stability problem for linear multistep methods," is associated with the A-stability concept and the second Dahlquist barrier.[14](https://doi.org/10.1007/bf01963532) David A. Pope's 1963 Communications of the ACM paper is an early exponential method for integrating ODEs.[15](https://doi.org/10.1145/366707.367592) G. Wanner, E. Hairer, and S. P. Nørsett's 1978 BIT Numerical Mathematics paper developed order stars, the theory clarifying how order and stability interact in implicit Runge–Kutta methods.[16](https://doi.org/10.1007/bf01932026) Marlis Hochbruck, Christian Lubich, and Hubert Selhofer's 1998 SIAM Journal on Scientific Computing paper treated [exponential integrators](https://www.edgepedia.org/exponential-integrator) for large systems,[17](https://doi.org/10.1137/s1064827595295337) and S. M. Cox and P. C. Matthews's 2002 Journal of Computational Physics paper introduced the ETDRK4 exponential time-differencing scheme.[18](https://doi.org/10.1006/jcph.2002.6995) J. C. Butcher's 1996 Applied Numerical Mathematics survey recounts the development of Runge–Kutta methods.[19](https://doi.org/10.1016/0168-9274%2895%2900108-5) In parallel-in-time integration, V. A. Dobrev and colleagues developed the two-level convergence theory for MGRIT in a 2017 SIAM Journal on Scientific Computing paper,[20](https://doi.org/10.1137/16m1074096) Alex C. Fish, Daniel R. Reynolds, and Steven B. Roberts presented implicit-explicit multirate infinitesimal stage-restart methods in a 2023 arXiv paper,[21](https://doi.org/10.48550/arxiv.2301.00865) and Martin J. Gander and Thibaut Lunet published the SIAM book *Time Parallel Time Integration* in 2024.[22](https://doi.org/10.1137/1.9781611978025)
## Variants
Geometric integrators preserve structure. The symplectic [Euler method](https://www.edgepedia.org/euler-method) updates momentum and position in a staggered way; the standard molecular-dynamics integrator is symplectic, and molecular simulators used symplectic methods for over twenty years before symplecticity entered numerical analysis.[1](https://unige.ch/~hairer/preprints/pcam-ode.pdf) For separable Hamiltonians \( H(q,p) = T(p) + V(q) \), Strang splitting reduces to the Störmer–Verlet method and Lie–Trotter splitting to the two symplectic Euler variants; compositions of symplectic maps stay symplectic, giving bounded energy error over long times.[23](https://www.cambridge.org/core/services/aop-cambridge-core/content/view/2C1BD434F8D2052E593A755351B3ACA3/S0962492923000077a.pdf/splitting_methods_for_differential_equations.pdf)
Exponential integrators solve the linear part exactly: \( y_{n+1} = y_n + h\, \varphi(h \cdot J_n) \cdot f_n \) is exact for \( f(y) = J \cdot y + c \).[1](https://unige.ch/~hairer/preprints/pcam-ode.pdf) IMEX schemes treat stiff terms implicitly and nonstiff terms explicitly; combining an explicit and an implicit Runge–Kutta method yields an IMEX ARK method, and in SDIRK methods each stage solve reuses the same matrix with a single diagonal \( \gamma \).[24](https://clima.github.io/ClimaTimeSteppers.jl/stable/algorithm_formulations/ode_solvers/) • [11](https://excalibur-neptune.github.io/Documents/_static/TN-03_AReviewTimeSteppingTechniquesPreconditioningHyperbolicAnisotropicEllipticProblem.pdf) Multirate methods such as MRIStep integrate fast and slow components with different steps.[12](https://drreynolds.github.io/files/atpesc-2022-tutorial.pdf)
## Applications
[Molecular dynamics](https://www.edgepedia.org/molecular-dynamics) runs on symplectic integrators whose long-time energy behavior is the point of the design.[1](https://unige.ch/~hairer/preprints/pcam-ode.pdf) Stiff chemistry calculations are the problem class that produced BDF methods.[5](http://www.scholarpedia.org/article/Linear_multistep_method) In fluid dynamics, direct numerical simulation of incompressible Navier–[Stokes flow](https://www.edgepedia.org/stokes-flow) past a circular cylinder has been demonstrated with a high-order stiffly stable splitting scheme inside a general time-stepping framework.[25](https://users.cs.utah.edu/~kirby/Publications/Kirby-53.pdf) Climate simulation uses dedicated stacks such as ClimaTimeSteppers, where IMEX Euler splits tendencies into explicit and implicit parts as a compromise between cost and stability.[24](https://clima.github.io/ClimaTimeSteppers.jl/stable/algorithm_formulations/ode_solvers/) In distributed-memory HPC, explicit solvers are usually preferred because they avoid the global communication of implicit solves.[26](https://openresearchsoftware.metajnl.com/articles/10.5334/jors.598)
## Limitations and alternatives
Known failure modes include order reduction: implicit Runge–Kutta methods applied to stiff problems can drop to observed order equal to the stage order or stage order plus 1.[1](https://unige.ch/~hairer/preprints/pcam-ode.pdf) Multistep methods carry a memory cost relative to Runge–Kutta methods, a concern at exascale, and no implicit multistep method is both A-stable and of order above two.[11](https://excalibur-neptune.github.io/Documents/_static/TN-03_AReviewTimeSteppingTechniquesPreconditioningHyperbolicAnisotropicEllipticProblem.pdf) IMEX splitting is not always clean: the split between regimes is not always clear, schemes may leak stiffness, and tight tolerances may force fully implicit treatment.[11](https://excalibur-neptune.github.io/Documents/_static/TN-03_AReviewTimeSteppingTechniquesPreconditioningHyperbolicAnisotropicEllipticProblem.pdf) Space-time methods, which discretize space and time simultaneously to overcome the limits of first discretizing in space, are an emerging and rapidly growing alternative.[27](https://publications.mfo.de/bitstream/handle/mfo/3934/OWR_2022_06.pdf?isAllowed=y&sequence=4) Parallel-in-time methods form another alternative family, grouped into shooting-type, waveform relaxation, space-time multigrid, and direct time-parallel approaches, with MGRIT reporting speedups up to 50× for linear parabolic problems while wrapping an existing sequential stepper.[28](https://www.cambridge.org/core/journals/acta-numerica/article/time-parallelization-for-hyperbolic-and-parabolic-problems/E580B25A44766E6729A65D7AF05B1198) • [29](https://www.osti.gov/servlets/purl/1236132)
## References
1. [Numerical solution of ordinary differential equations (PCAM chapter, Hairer & Lubich)](https://unige.ch/~hairer/preprints/pcam-ode.pdf)
2. [Numerical Methods for solving ODEs and PDEs (EMSE course notes)](https://www.emse.fr/~bonnefoy/Public/MetNum-EMSE.pdf)
3. [Some simple time-stepping methods (Schöberl, TU Wien)](https://jschoeberl.github.io/IntroSC/ODEs/firstmethods.html)
4. [DuMux: Benchmark Time-Stepping Methods](https://dumux.org/docs/doxygen/master/benchmark-timestepping-methods.html)
5. [Linear multistep method (Scholarpedia, Hairer/Lubich/Wanner)](http://www.scholarpedia.org/article/Linear_multistep_method)
6. [Runge–Kutta Implementation, Adaptive and Multirate Time Integration (Reynolds, 2025 tutorial)](https://drreynolds.github.io/files/viasm-2025-tutorial.pdf)
7. [Review of Time-Stepping Algorithms for ODEs (Brown University course notes)](https://bpb-us-w2.wpmucdn.com/sites.brown.edu/dist/1/376/files/2022/04/Chapter-3-Review-of-Time-Stepping-Algorithms-for-ODEs.pdf)
8. [Finite Difference Methods for Differential Equations (ODEs and RK schemes)](https://djps.github.io/docs/numericalmethods/part3/odes/)
9. [Numerical Methods for the Solution of Differential Equations (University of Vienna)](https://www.csc.univie.ac.at/files/Numerical_methods_for_differential_equations_W2015/main_1.pdf)
10. [Efficient solvers for time-dependent problems: IMEX, LATIN, PARAEXP, PARAREAL (Springer review)](https://link.springer.com/article/10.1186/s40323-016-0063-y)
11. [A Review of Time Stepping Techniques and Preconditioning for Hyperbolic Problems](https://excalibur-neptune.github.io/Documents/_static/TN-03_AReviewTimeSteppingTechniquesPreconditioningHyperbolicAnisotropicEllipticProblem.pdf)
12. [Time Integration with SUNDIALS (ATPESC 2022 tutorial)](https://drreynolds.github.io/files/atpesc-2022-tutorial.pdf)
13. [C. F. Curtiss, J. O. Hirschfelder (1952). Integration of Stiff Equations. Proceedings of the National Academy of Sciences.](https://doi.org/10.1073/pnas.38.3.235)
14. [Germund G. Dahlquist (1963). A special stability problem for linear multistep methods. BIT Numerical Mathematics.](https://doi.org/10.1007/bf01963532)
15. [David A. Pope (1963). An exponential method of numerical integration of ordinary differential equations. Communications of the ACM.](https://doi.org/10.1145/366707.367592)
16. [G. Wanner, E. Hairer, S. P. Nørsett (1978). Order stars and stability theorems. BIT Numerical Mathematics.](https://doi.org/10.1007/bf01932026)
17. [Marlis Hochbruck, Christian Lubich, Hubert Selhofer (1998). Exponential Integrators for Large Systems of Differential Equations. SIAM Journal on Scientific Computing.](https://doi.org/10.1137/s1064827595295337)
18. [S.M. Cox, P.C. Matthews (2002). Exponential Time Differencing for Stiff Systems. Journal of Computational Physics.](https://doi.org/10.1006/jcph.2002.6995)
19. [A history of Runge-Kutta methods (Applied Numerical Mathematics, 1996)](https://doi.org/10.1016/0168-9274%2895%2900108-5)
20. [V. A. Dobrev and colleagues (2017). Two-Level Convergence Theory for Multigrid Reduction in Time (MGRIT). SIAM Journal on Scientific Computing.](https://doi.org/10.1137/16m1074096)
21. [Fish, Alex C., Reynolds, Daniel R., Roberts, Steven B. (2023). Implicit-Explicit Multirate Infinitesimal Stage-Restart Methods. arXiv (Cornell University).](https://doi.org/10.48550/arxiv.2301.00865)
22. [Martin J. Gander, Thibaut Lunet (2024). Time Parallel Time Integration. Society for Industrial and Applied Mathematics eBooks.](https://doi.org/10.1137/1.9781611978025)
23. [Splitting methods for differential equations (Acta Numerica review)](https://www.cambridge.org/core/services/aop-cambridge-core/content/view/2C1BD434F8D2052E593A755351B3ACA3/S0962492923000077a.pdf/splitting_methods_for_differential_equations.pdf)
24. [ODE Solvers · ClimaTimeSteppers.jl documentation](https://clima.github.io/ClimaTimeSteppers.jl/stable/algorithm_formulations/ode_solvers/)
25. [A Generic Framework for Time-Stepping PDEs (Kirby et al.)](https://users.cs.utah.edu/~kirby/Publications/Kirby-53.pdf)
26. [Integrating Odeint Time Stepping into OpenFPM (JOSS/JORS)](https://openresearchsoftware.metajnl.com/articles/10.5334/jors.598)
27. [Space-Time Methods for Time-Dependent PDEs (Oberwolfach Workshop Report, 2022)](https://publications.mfo.de/bitstream/handle/mfo/3934/OWR_2022_06.pdf?isAllowed=y&sequence=4)
28. [Time parallelization for hyperbolic and parabolic problems (Acta Numerica)](https://www.cambridge.org/core/journals/acta-numerica/article/time-parallelization-for-hyperbolic-and-parabolic-problems/E580B25A44766E6729A65D7AF05B1198)
29. [Multigrid Reduction in Time for Nonlinear Parabolic Problems (OSTI report)](https://www.osti.gov/servlets/purl/1236132)
---
*Topic: Encyclopedia › Physical world and mathematics › Mathematics and statistics › Analysis and mathematical models › Numerical analysis and computation › Time integration methods*
*Initially written Sep 29, 2026 · Edited: Oct 15, 2026*
*Copyright 2026 EdgeChat AI, a subsidiary of Biostate AI.*
License: Edgepedia Community License 1.0, https://www.edgepedia.org/edgepedia/license