The *fixed-leading coefficient* form of BDF, the choice in all SUNDIALS integrators can be derived (and implemented) most conveniently via interpolating polynomials. The derivation for BDFq can be obtained in the following 10 steps.

Suppose we have computed $y_j \approx y(t_j), j = 1,2,....,n-1$. Let the past step sizes be $h_j = t_j - t_{j-1}, j=1,2,...,n-1$ and $h_n = t_n - t_{n-1}$ be the current step size. We wish to compute $y_n \approx y(t_n)$.

**Step 1**

Build a predictor polynomial $\omega_p$ of order $q$ that satisfies the following $(q+1)$ conditions:

In other words, $\omega_p$ interpolates the past $q$ computed values plus the derivative at the last computed step.

**Step 2**

Build a corrector polynomial $\omega_c$ of order $q$ that satisfies the following $(q+1)$ conditions:

In other words, $\omega_c$ passes through the (unknown at this point) solution at $t_n$ and coincides with the predictor polynomial at $q$ past equidistant points (separated by the current step size $h_n$).

**Step 3**

It is easy to see that the $q$-th order FLC BDF formula can be written simply as

**Step 4**

Since both $\omega_p$ and $\omega_c$ are $q$-th order polynomial, let's relate them through a new polynomial $c(t)$ (which we will define later) as

It's obvious that, once plugged into the ODE, the above BDF formula will lead to a nonlinear equation in the desired solution $y_n$. We will use some form of iteration to solve this nonlinear equation for $y$ at $t_n$ (recall that this is, by construction, $\omega_c(t_n)$ ) and we will use the predicted value at $t_n$ as our initial guess.

**Step 5**

Denote then

Note that $y_{n(m)}, m=1.2,...$ will be successive iterates of the nonlinear solver (and $y_n$ will be $y_{n(M)}$ where $M$ is the iteration at which we declare convergence).

**Step 6**

With the above definitions we have

**Step 7**

Now use the relations (1) and (2) to derive conditions on $c(t)$. In other words, evaluate (6) at $t_n - i h_n, i=0,1,...,q$ and use (1) and (2) to get:

which are $(q+1)$ conditions uniquely defining the $q$-th order polynomial $c(t)$.

**Step 8**

Using the relation (6), the BDF formula in (3) can now be written as:

or, using (5)

(9)All we need to do now to complete the derivation of the BDF formula is find the derivative of $c(t)$ at $t_n$.

**Step 9**

Given the conditions (7), i.e. the roots of the $q$-th order polynomial $c(t)$, we have:

And at $t=t_n$, we get

(11)Therefore $K = {1}/{q!h_n^q}$ and $c(t)$ is completely determined.

**Step 10**

Differentiate (10) and evaluate at $t=t_n$ to get

or

(13)where we define $\beta_0 = 1/(1+1/2+...+1/q)$.

All done. Multiply (9) by $h_n$ and use (13) to get the BDF formula as:

(14)or

(15)or

(16)