CVODE, CVODES, and IDA all use Backward Differentiation Formulas methods (CVODE and CVODES also provide the option of using AdamsMoulton methods for the solution of nonstiff problems).
Table of Contents

Basic formulas
Denote by $t_n$ and $y_n = y(t_n)$ the discrete values of the independent variable $t$ and solution $y(t)$ and let $h = t_n  t_{n1}$ be the stepsize.
BDF is one of many Linear Multistep Methods (LMM). It is called BDF because it differentiates the solution $y$ using past (backward) values; i.e.,
(1)The simplest BDF method (BDF1) is Backward Euler (BE):
(2)The fixedstep BDF of order $q$ (BDFq) is
(3)where the coefficients $\alpha_i$ and $\beta_0$ depend only on the method order $q$.
For an ODE in explicit form; i.e. $\dot y = f(t,y)$, the BDFq method gives $y_n$ as the solution of
(4)For an ODE in implicit form or a DAE; i.e. $F(t,y,\dot y) = 0$, the BDFq method gives $y_n$ as the solution of
(5)Symbolic operator derivation
For any sequence $\{x_n\}$, define the operators
(6)It is easy to develop an algebra of operators by observing that the following relationships hold:
(7)where the inverse increment operator is defined by
(8)A given sequence may or may not be the set of discrete values of a smooth functino of $t$. When it is, we can define a differentiation operator, $Df_n = \dot f(t_n)$ and write the infinite Taylor series expansion of $f(t)$ at $t = t_n$ as
(9)In other words, the increment and differentiation operators are related through
(10)Next, using (7), we get $E = (1\nabla)^{1}$ and therefore
(11)which, applied to $y(t)$ gives:
(12)to obtain the BDFq method, simply truncate the series (12) at the $\nabla^q$ term. Through simple manipulations, it is easy to obtain
(13)Varying the stepsize
There are three possible ways of allowing a variable stepsize: interpolated fixedstep BDF, fully variablestep BDF, and fixedleading coefficient BDF.
(1) In the interpolated fixedstep approach, one uses the fixedstep BDF method with stepsize $h$, but interpolates at a spacing equal to the new stepsize $h'$. However, the behavior of the resulting method can be unstable if the stepsize changes too frequently.
(2) A fully variablestep BDF can be obtained, by allowing the method coefficients to also depend on the stepsize $h$, as
(14)Unlike the method (3), the coefficients $\alpha_{n,i}$ and $\beta_{n,0}$ depend on the method order $q$ and the last $q$ stepsizes $h_n , h_{n1} , \ldots , h_{nq+1}$. Denoting $h = \max \{ h_n , h_{n1} , \ldots , h_{nq+1} \}$, the method coefficients are uniquely determined by
(15)The disadvantage of this approach is that the variation of the quantity $h_n \beta_{n,0}$ makes it harder to keep the Newton iteration matrix (for the solution of the resulting nonlinear system) current.
(3) The fixedleading coefficient form of BDF is a compromise between the above two choices in that it has the fixedstep value (13) of $\beta_{n,0} \equiv \beta_0 = 1 / \sum_{j=1}^q j^{1}$ but at the cost of introducing an additional term $\beta_{n,1}} \dot y_{n1}$:
(16)As expected, the coefficient $\beta_{n,1}$ vanishes when the stepsizes are equal. The fixedleading coefficient form of BDF is most conveniently derived via interpolating polynomials. This derivation is explained in details in a separate article.
The main advantages of FLC BDF is that it does not suffer from the unstable behavior of the interpolated fixedstep method and, since the only variation in the quantity $h_n \beta_0$ is now through the stepsize $h_n$ only, the Newton iteration matrix can be reused for more steps than in a fully variablestep approach.
The fixedleading coefficient form of BDF is the choice in all SUNDIALS integrators (CVODE, CVODES, IDA).
Example: BDF2
Fixedstep
(17)Variablestep
(18)Fixedleading coefficient
(19)