Numerical Methods for ODE I have been working through Iserles books A first course in Numerical Analysis of Ordinary Differential Equations. 
I am trying to figure out how to prove that the error constant for an s-step BDF is -beta/(s+1) 
I have tried using the characteristic polynomials to no avail. 
does anyone have any ideas? 
 A: Multistep methods in general follow the scheme
$$
\sum a_k y(x-kh)\approx h\sum b_k f(x-kh, y(x-kh))
$$
For an exact solution we have that $f(x,y(x))=y'(x)$ and the order condition can be then written as
$$
\sum a_k y(x-kh)=h\sum b_k y'(x-kh) + O(h^{p+1})
$$
The coefficients inside this formula have to be valid for all differentiable functions. It is sufficient to explore it for all analytical functions. Then one can write the shift operator for $y(x-h)$ as $e^{-hD}$ with $D=\frac{d}{dx}$, so that the order condition translates to
$$
\sum a_k e^{-khD} = hD\sum b_k e^{-khD} + O(h^{p+1})
$$
Now write $e^{-hD}=1-Z$, $hD=-\ln(1-Z)$ and we get an equation in truncated power series
$$
\sum a_k(1-Z)^k=-\ln(1-Z)\sum b_k(1-Z)^k + O(Z^{p+1})
$$

For BDF methods one uses (for the moment) $b_0=1$ and $b_k=0$ for $k>1$ so that
$$
a(1-Z)=Z+\frac12Z^2+…+\frac1pZ^p
$$
and the error terms is $\frac1{p+1}Z^{p+1}+O(Z^{p+2})$. Since $Z=1-e^{-hD}=hD+O(h^2)$ one gets $Z^{p+1}y=h^{p+1}y^{(p+1)}(x)+O(h^{p+2})$. In the coefficients $a_k$ we get
$$
\sum_{k=0}^p a_kq^k = (1-q)+\frac12(1-q)^2+…+\frac1p(1-q)^p\\
=\sum_{k=0}^p q^k\sum_{j=k}^p\binom{j}{k}\frac1j
$$
Obviously, the constant coefficient is $a_0=1+\frac12+…+\frac1p$ and by multiplying all coefficients with $β=(1+\frac12+…+\frac1p)^{-1}$ one gets a method with $a_0=1$, $b_0=β$ and leading error term 
$$
\frac{β}{p+1}h^{p+1}y^{(p+1)}(x)+O(h^{p+2}).
$$
