  In this paper, we develop a numerical scheme for the space-time fractional
parabolic equation, i.e., an equation involving a fractional time derivative
and a fractional spatial operator. Both the initial value problem and the
non-homogeneous forcing problem (with zero initial data) are considered. The
solution operator $E(t)$ for the initial value problem can be written as a
Dunford-Taylor integral involving the Mittag-Leffler function $e_{\alpha,1}$
and the resolvent of the underlying (non-fractional) spatial operator over an
appropriate integration path in the complex plane. Here $\alpha$ denotes the
order of the fractional time derivative. The solution for the non-homogeneous
problem can be written as a convolution involving an operator $W(t)$ and the
forcing function $F(t)$.
  We develop and analyze semi-discrete methods based on finite element
approximation to the underlying (non-fractional) spatial operator in terms of
analogous Dunford-Taylor integrals applied to the discrete operator. The space
error is of optimal order up to a logarithm of $1/h$. The fully discrete method
for the initial value problem is developed from the semi-discrete approximation
by applying an exponentially convergent sinc quadrature technique to
approximate the Dunford-Taylor integral of the discrete operator and is free of
any time stepping.
  To approximate the convolution appearing in the semi-discrete approximation
to the non-homogeneous problem, we apply a pseudo midpoint quadrature. This
involves the average of $W_h(s)$, (the semi-discrete approximation to $W(s)$)
over the quadrature interval. This average can also be written as a
Dunford-Taylor integral. We first analyze the error between this quadrature and
the semi-discrete approximation. To develop a fully discrete method, we then
introduce sinc quadrature approximations to the Dunford-Taylor integrals for
computing the averages.
