We describe an algorithm for solving systems of parabolic partial differential equations in three space dimensions.  Approximate solutions are determined using Galerkin's method with a high-order, piecewise polynomial hierarchical basis in space and the diferential algebraic equations code DASPK in time.  On three test problems, a variant of DASPK that computes Jacobian vector products with a stored analytic Jacobian is compared with the default matrix-free approximation of the Jacobian vector product.  Several preconditioners are compared including incomplete Cholesky, banded LU, and element-by-element factorizations.