Solution (source code)

= Solution

Interpolate through the $s+1$ values at $t_n,\ldots,t_{n+s}$, and differentiate the interpolation <polynomial> at its newest node. In the backward coordinate $\tau=(t-t_{n+s})/h$, the Newton form is
$$
P(\tau)=\sum_{j=0}^s\frac{\tau(\tau+1)\cdots(\tau+j-1)}{j!}\,\nabla^j y_{n+s},
$$
with the $j=0$ factor interpreted as one and $\nabla y_m=y_m-y_{m-1}$. For $j\ge1$, the derivative of the factor at zero is $(j-1)!/j!=1/j$. Thus the <backward differentiation formula> is
$$
\sum_{j=1}^s\frac1j\nabla^j y_{n+s}=h f_{n+s}.
$$
If $w$ denotes the forward-shift variable, $\nabla^j y_{n+s}$ corresponds to $w^{s-j}(w-1)^j$. Consequently the requested <polynomial> is
$$
\boxed{\rho(w)=\sum_{j=1}^s\frac1j w^{s-j}(w-1)^j
=w^s\sum_{j=1}^s\frac{(1-w^{-1})^j}{j}.}
$$
The first form is a <polynomial> valid also at $w=0$. The second form is convenient for relating it to the truncated series of $\log w$. Differentiation of the degree-$s$ interpolant has local derivative error $O(h^s)$, giving order $s$ for smooth solutions.