  We propose and study a class of numerical schemes to approximate time
fractional differential equations. The methods are based on the approximation
of the Caputo fractional derivative by continuous piecewise polynomials, which
is strongly related to the backward differentiation formulae for the
integer-order case. We investigate their theoretical properties, such as the
local truncation error and global error analyses with respect to a sufficiently
smooth solution, and the numerical stability in terms of the stability region
and $A(\frac{\pi}{2})$-stability by refining the technique proposed in
\cite{LubichC:1986b}. Numerical experiments are given to verify the theoretical
investigations.
