  This article extends the semidiscrete maximal $L^p$-regularity results in
[27] to multistep fully discrete finite element methods for parabolic equations
with more general diffusion coefficients in $W^{1,d+\beta}$, where $d$ is the
dimension of space and $\beta>0$. The maximal angles of $R$-boundedness are
characterized for the analytic semigroup $e^{zA_h}$ and the resolvent operator
$z(z-A_h)^{-1}$, respectively, associated to an elliptic finite element
operator $A_h$. Maximal $L^p$-regularity, optimal $\ell^p(L^q)$ error estimate,
and $\ell^p(W^{1,q})$ estimate are established for fully discrete finite
element methods with multistep backward differentiation formula.
