# Backward differentiation formula

The backward differentiation formula (BDF) is a family of implicit linear multistep methods for computing numerical solutions of ordinary differential equations, designed for stiff systems.<sup>[1](http://www.scholarpedia.org/article/Backward_differentiation_formulas)</sup> Each step replaces the time derivative by a backward difference through recent solution values, leaving one implicit equation to solve per step; this implicitness is essential for arbitrarily stiff systems.<sup>[1](http://www.scholarpedia.org/article/Backward_differentiation_formulas)</sup> The first use of BDF methods is credited to Curtiss and Hirschfelder's 1952 paper on stiff equations,<sup>[2](https://doi.org/10.1073/pnas.38.3.235)</sup> and the family now underlies standard stiff solvers in scientific computing, including CVODE, DASSL, and MATLAB's ode15s.<sup>[3](https://sundials.readthedocs.io/en/v7.1.0/cvode/Mathematics_link.html)</sup><sup> • </sup><sup>[4](https://www.osti.gov/servlets/purl/5882821)</sup><sup> • </sup><sup>[5](http://bicycle.tudelft.nl/schwab/TAM674/SR97.pdf)</sup>

| Key fact | Detail |
| --- | --- |
| Method class | Implicit linear multistep methods of orders 1 through 6, for stiff ODEs and some differential-algebraic equations<sup>[1](http://www.scholarpedia.org/article/Backward_differentiation_formulas)</sup> |
| Defining formula | \( h \cdot y'_{n} = \sum_{j=0}^{k} \alpha_{k,j} \, y_{n-j} \) at fixed step \( h \)<sup>[1](http://www.scholarpedia.org/article/Backward_differentiation_formulas)</sup> |
| Coefficients | Taken from the Lagrange interpolation polynomial through recent values; BDF1 is the backward Euler method<sup>[6](https://sofa-framework.github.io/doc/components/odesolver/backward/bdfodesolver/)</sup> |
| Stability | A-stable at orders 1 and 2; A(α)-stable at orders 3 to 6, with the angle shrinking from 90° to about 17.84°<sup>[3](https://sundials.readthedocs.io/en/v7.1.0/cvode/Mathematics_link.html)</sup><sup> • </sup><sup>[7](https://www.cse.uoi.gr/~akrivis/AK2a.pdf)</sup> |
| Order barrier | For \( k \geq 7 \) the schemes are not even 0-stable, so no BDF beyond order 6 is usable<sup>[8](https://www.opuscula.agh.edu.pl/vol26/2/art/opuscula_math_2615.pdf)</sup> |
| Stiff stability | Gear's stiffness limits \( \delta \): 0 for orders 1 and 2, then −0.1, −0.7, −2.4, −6.1 for orders 3 to 6<sup>[9](https://arxiv.org/html/1101.4434)</sup> |
| Typical software | CVODE (orders 1 to 5, fixed-leading coefficient form), DASSL, ode15s<sup>[3](https://sundials.readthedocs.io/en/v7.1.0/cvode/Mathematics_link.html)</sup><sup> • </sup><sup>[4](https://www.osti.gov/servlets/purl/5882821)</sup><sup> • </sup><sup>[5](http://bicycle.tudelft.nl/schwab/TAM674/SR97.pdf)</sup> |

## How it works

At a constant step size \( h \), the order-\( k \) BDF for a step to \( (t_{n+1}, y_{n+1}) \) is

\[ h \cdot y'_{n+1} = \sum_{j=0}^{k} \alpha_{k,j} \, y_{n+1-j} \]

with coefficients taken from the Lagrange interpolation polynomial through the recent solution values; the order equals the number of previous time steps the polynomial requires, and order 1 is the backward [Euler method](https://www.edgepedia.org/euler-method).<sup>[1](http://www.scholarpedia.org/article/Backward_differentiation_formulas)</sup><sup> • </sup><sup>[6](https://sofa-framework.github.io/doc/components/odesolver/backward/bdfodesolver/)</sup> Because \( y_{n} \) appears on both sides of the equation the method is implicit, and it uses only one implicit evaluation of the right-hand side of the ODE per step.<sup>[1](http://www.scholarpedia.org/article/Backward_differentiation_formulas)</sup><sup> • </sup><sup>[8](https://www.opuscula.agh.edu.pl/vol26/2/art/opuscula_math_2615.pdf)</sup>

**Stability is the design target.** Orders 1 and 2 are A-stable: for any complex \( \lambda \) in the open left half-plane, the method is unconditionally stable at any step size for the model problem \( y' = \lambda y \).<sup>[3](https://sundials.readthedocs.io/en/v7.1.0/cvode/Mathematics_link.html)</sup> Orders 3 through 6 are A(α)-stable with maximum angles \( \vartheta_{1} = \vartheta_{2} = 90^{\circ} \), \( \vartheta_{3} \approx 86.03^{\circ} \), \( \vartheta_{4} \approx 73.35^{\circ} \), \( \vartheta_{5} \approx 51.84^{\circ} \), and \( \vartheta_{6} \approx 17.84^{\circ} \).<sup>[7](https://www.cse.uoi.gr/~akrivis/AK2a.pdf)</sup> The same orders are stiffly stable in Gear's sense, with stiffness limits \( \delta = 0 \) for orders 1 and 2 and \( \delta = -0.1, -0.7, -2.4, -6.1 \) for orders 3 to 6.<sup>[9](https://arxiv.org/html/1101.4434)</sup> For \( k \geq 7 \) the schemes are not even 0-stable,<sup>[8](https://www.opuscula.agh.edu.pl/vol26/2/art/opuscula_math_2615.pdf)</sup> and the order-7 stability contour crosses the negative real axis, making order 7 and higher of no value; order 6, though still A(α)-stable, has so tight an angle that it is rarely implemented.<sup>[1](http://www.scholarpedia.org/article/Backward_differentiation_formulas)</sup><sup> • </sup><sup>[8](https://www.opuscula.agh.edu.pl/vol26/2/art/opuscula_math_2615.pdf)</sup>

## How it is done

Quality codes start BDF integration with Adams methods, which are more effective on the initial transients common to stiff problems and whose order-selection mechanism reduces order appropriately on problems with eigenvalues near the imaginary axis.<sup>[10](https://stewart.sdsu.edu/mypubs/StewartJCAM_Avoiding_1-s2.0-037704279090017T-main.pdf)</sup> The implicit equation at each step is then solved by a Newton or modified Newton iteration requiring the solution of linear systems.<sup>[4](https://www.osti.gov/servlets/purl/5882821)</sup><sup> • </sup><sup>[3](https://sundials.readthedocs.io/en/v7.1.0/cvode/Mathematics_link.html)</sup>

**Error and order control.** CVODE controls local error through the asymptotic estimate \( \mathrm{LTE} = C \cdot h^{q+1} \cdot y^{(q+1)} + O(h^{q+2}) \) at order \( q \) and step size \( h \), redoing steps with a reduced step size when the error test fails; order is varied dynamically between 1 and 5, with changes considered only after \( q+1 \) steps at order \( q \), and then only to \( q-1 \) or \( q+1 \).<sup>[3](https://sundials.readthedocs.io/en/v7.1.0/cvode/Mathematics_link.html)</sup> For a fixed step size \( h \), a \( k \)-step BDF with \( k < 7 \) is convergent of order \( k \) if all initial values are accurate to \( O(h^{k}) \) and the Newton iteration converges to \( O(h^{k+1}) \); this extends to variable step sizes provided the implementation is stable for ODEs.<sup>[11](https://www.cs.usask.ca/~spiteri/M314/notes/AP/chap10.pdf)</sup> CVODE additionally includes the STALD (stability limit detection) algorithm to protect against potentially unstable BDF behavior.<sup>[3](https://sundials.readthedocs.io/en/v7.1.0/cvode/Mathematics_link.html)</sup>

## Origin

The first use of BDF methods is credited to C. F. Curtiss and J. O. Hirschfelder's paper "Integration of Stiff Equations", published in the [Proceedings of the National Academy of Sciences](https://www.edgepedia.org/proceedings-of-the-national-academy-of-sciences) in 1952 (38(3):235–243).<sup>[2](https://doi.org/10.1073/pnas.38.3.235)</sup> The paper addresses equations arising in chemical kinetics, electrical circuit theory, and missile guidance that are exceedingly difficult to solve by ordinary numerical procedures, and proposes a forward interpolation (backward differentiation) scheme; the same source is credited with coining the notion of stiff differential equations, and one of the first definitions of stiffness: equations on which certain implicit methods outperform explicit ones such as Euler or Adams methods.<sup>[12](https://unige.ch/~hairer/preprints/pcam-ode.pdf)</sup><sup> • </sup><sup>[9](https://arxiv.org/html/1101.4434)</sup> Henrici (1962) discussed methods based on differentiation but dismissed them as less accurate than the corresponding Adams-Moulton formulas for non-stiff equations.<sup>[1](http://www.scholarpedia.org/article/Backward_differentiation_formulas)</sup>

C. W. Gear's [Algorithm](https://www.edgepedia.org/algorithm) 407, DIFSUB, published in Communications of the ACM in 1971,<sup>[13](https://doi.org/10.1145/362566.362573)</sup> was the first widely used code for stiff problems and formalized the method often called Gear's method.<sup>[12](https://unige.ch/~hairer/preprints/pcam-ode.pdf)</sup> The DASSL report extended BDF to differential-algebraic systems, approximating the derivative with BDF of order one to five and choosing order and step size from solution behavior.<sup>[4](https://www.osti.gov/servlets/purl/5882821)</sup>

## Variants

**Step-size formulations.** Variable-step codes represent the method's history in modified divided differences or the Nordsieck array.<sup>[1](http://www.scholarpedia.org/article/Backward_differentiation_formulas)</sup> DASSL uses the fixed-leading coefficient form (after Jackson and Sacks-Davis), which tends to be more stable than the fixed coefficient formulas used in LSODI and more efficient than the variable coefficient formulas used in VODE.<sup>[4](https://www.osti.gov/servlets/purl/5882821)</sup> CVODE likewise uses variable-order, variable-step BDF in fixed-leading coefficient form of orders 1 to 5, with coefficients determined by the method type, its order, the recent step-size history, and the normalization \( \alpha_{n,0} = -1 \).<sup>[3](https://sundials.readthedocs.io/en/v7.1.0/cvode/Mathematics_link.html)</sup>

**DAE and split forms.** BDF can be applied directly to some differential-algebraic equations of the form \( F(y', y, t) = 0 \),<sup>[1](http://www.scholarpedia.org/article/Backward_differentiation_formulas)</sup> and most BDF-based software packages solve fully implicit index-1 DAEs.<sup>[11](https://www.cs.usask.ca/~spiteri/M314/notes/AP/chap10.pdf)</sup> IMEX (semi-implicit) BDF schemes apply the implicit BDF to the stiff part of a split system and treat the non-stiff part explicitly, with accuracy up to fourth order; the OrdinaryDiffEq.jl implementations include SBDF2 (recommended), SBDF3, SBDF4, and an experimental adaptive-order variant.<sup>[14](https://github.com/SciML/OrdinaryDiffEq.jl/blob/master/docs/src/imex/IMEXBDF.md)</sup> K. Stewart's 1990 paper in the Journal of Computational and Applied Mathematics modeled the stability of these semi-implicit formulas, treating a prediction followed by a fixed number of Newton corrections made with an inexact Jacobian matrix.<sup>[15](https://doi.org/10.1016/0377-0427%2890%2990050-a)</sup>

**Software.** MATLAB's ODE suite implements the BDF-based stiff solver ode15s.<sup>[5](http://bicycle.tudelft.nl/schwab/TAM674/SR97.pdf)</sup> TensorFlow Probability's JAX substrate provides a variable-step, variable-order BDF integrator with order between 1 and 5, following Shampine and Reichelt (1997), with gradients computed by the adjoint sensitivity method.<sup>[16](https://www.tensorflow.org/probability/api_docs/python/tfp/substrates/jax/math/ode/BDF)</sup> VODE is a variable-step BDF method in fixed-leading coefficient form, and CVODE is a C-language successor of VODE within the SUNDIALS suite;<sup>[21](https://sundials.readthedocs.io/en/v7.5.0/cvode/Organization_link.html)</sup> multistep BDF methods solve a single implicit equation per step, whereas Rosenbrock and (E)SDIRK methods must solve a few.<sup>[17](https://www.stochasticlifestyle.com/differences-between-methods-for-solving-stiff-odes/)</sup>

## Applications

DifferentialEquations.jl documents BDF as the preferred choice for very large stiff systems with more than 1000 equations, where other implicit methods become computationally expensive, and lists reaction-diffusion systems, chemical kinetics, circuit simulation, and parabolic PDEs after spatial discretization as recommended use cases.<sup>[18](https://docs.sciml.ai/DiffEqDocs/stable/api/ordinarydiffeq/implicit/BDF/)</sup> For large stiff systems where direct linear solvers are infeasible, combining a BDF integrator with a preconditioned Krylov method is described as a powerful tool.<sup>[3](https://sundials.readthedocs.io/en/v7.1.0/cvode/Mathematics_link.html)</sup>

## Limitations and alternatives

Orders 3 through 6 have an exit angle below 60° when eigenvalues lie near the imaginary axis, which causes instability there;<sup>[10](https://stewart.sdsu.edu/mypubs/StewartJCAM_Avoiding_1-s2.0-037704279090017T-main.pdf)</sup> the CVODE documentation likewise notes that at orders 3 to 5 the methods are not A-stable, with a region of instability near the imaginary axis that grows with order.<sup>[3](https://sundials.readthedocs.io/en/v7.1.0/cvode/Mathematics_link.html)</sup> Because only orders 1 and 2 are L-stable, BDF loses effectiveness at high orders on more than semi-stiff problems and must shrink the time step.<sup>[18](https://docs.sciml.ai/DiffEqDocs/stable/api/ordinarydiffeq/implicit/BDF/)</sup> Radau IIA methods, implicit Runge-Kutta schemes that are both A-stable and L-stable, are a leading alternative; Hairer and Wanner's RADAU implementation with a variable order strategy is presented as an alternative to BDF codes.<sup>[19](https://www.sciencedirect.com/science/article/pii/S037704279900134X)</sup><sup> • </sup><sup>[20](https://doi.org/10.1016/s0377-0427%2899%2900134-x)</sup>

## References

1. [Backward differentiation formulas - Scholarpedia](http://www.scholarpedia.org/article/Backward_differentiation_formulas)
2. [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)
3. [CVODE Mathematical Considerations (SUNDIALS documentation)](https://sundials.readthedocs.io/en/v7.1.0/cvode/Mathematics_link.html)
4. [Description of DASSL: A Differential/Algebraic System Solver (Petzold)](https://www.osti.gov/servlets/purl/5882821)
5. [The MATLAB ODE Suite (Shampine & Reichelt)](http://bicycle.tudelft.nl/schwab/TAM674/SR97.pdf)
6. [BDFOdeSolver - SOFA Documentation](https://sofa-framework.github.io/doc/components/odesolver/backward/bdfodesolver/)
7. [Maximum angles of A(ϑ)-stability of backward difference formulae (Akrivis & Katsoprinakis)](https://www.cse.uoi.gr/~akrivis/AK2a.pdf)
8. [A uniform quantitative stiff stability estimate for BDF schemes (Opuscula Mathematica)](https://www.opuscula.agh.edu.pl/vol26/2/art/opuscula_math_2615.pdf)
9. [Analyzing the absolute stability region of implicit methods of solving ODEs](https://arxiv.org/html/1101.4434)
10. [Avoiding instability in BDF methods (Stewart, J. Comput. Appl. Math.)](https://stewart.sdsu.edu/mypubs/StewartJCAM_Avoiding_1-s2.0-037704279090017T-main.pdf)
11. [Chapter 10: BDF and Multistep Methods (course notes, U. Saskatchewan)](https://www.cs.usask.ca/~spiteri/M314/notes/AP/chap10.pdf)
12. [Hairer, Numerical analysis lecture notes (pcam-ode), Section 4.1 BDF methods](https://unige.ch/~hairer/preprints/pcam-ode.pdf)
13. [C. W. Gear (1971). Algorithm 407: DIFSUB for solution of ordinary differential equations [D2]. Communications of the ACM.](https://doi.org/10.1145/362566.362573)
14. [IMEX BDF methods (OrdinaryDiffEq.jl documentation)](https://github.com/SciML/OrdinaryDiffEq.jl/blob/master/docs/src/imex/IMEXBDF.md)
15. [A model for stability of the semi-implicit backward differentiation formulas (Journal of Computational and Applied Mathematics, 1990)](https://doi.org/10.1016/0377-0427%2890%2990050-a)
16. [tfp.substrates.jax.math.ode.BDF | TensorFlow Probability](https://www.tensorflow.org/probability/api_docs/python/tfp/substrates/jax/math/ode/BDF)
17. [Differences Between Methods for Solving Stiff ODEs - Stochastic Lifestyle (Christopher Rackauckas)](https://www.stochasticlifestyle.com/differences-between-methods-for-solving-stiff-odes/)
18. [OrdinaryDiffEqBDF · DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/api/ordinarydiffeq/implicit/BDF/)
19. [Stiff differential equations solved by Radau methods (Hairer & Wanner, J. Comput. Appl. Math.)](https://www.sciencedirect.com/science/article/pii/S037704279900134X)
20. [Stiff differential equations solved by Radau methods (Journal of Computational and Applied Mathematics, 1999)](https://doi.org/10.1016/s0377-0427%2899%2900134-x)
21. [Organization link (sundials.readthedocs.io)](https://sundials.readthedocs.io/en/v7.5.0/cvode/Organization_link.html)

---
*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 · Reviewed: Sep 30, 2026 · Edited: Sep 30, 2026 · Last review: Sep 30, 2026*

*Copyright 2026 EdgeChat AI, a subsidiary of Biostate AI.*

License: Edgepedia Community License 1.0, https://www.edgepedia.org/edgepedia/license