  We propose a piecewise-linear, time-stepping discontinuous Galerkin method to
solve numerically a time fractional diffusion equation involving Caputo
derivative of order $\mu\in (0,1)$ with variable coefficients. For the spatial
discretization, we apply the standard piecewise linear continuous Galerkin
method. Well-posedness of the fully discrete scheme and error analysis will be
shown. For a time interval~$(0,T)$ and a spatial domain~$\Omega$, our analysis
suggest that the error in $L^2\bigr((0,T),L^2(\Omega)\bigr)$-norm is of order
$O(k^{2-\frac{\mu}{2}}+h^2)$ (that is, short by order $\frac{\mu}{2}$ from
being optimal in time) where $k$ denotes the maximum time step, and $h$ is the
maximum diameter of the elements of the (quasi-uniform) spatial mesh. However,
our numerical experiments indicate optimal $O(k^{2}+h^2)$ error bound in the
stronger $L^\infty\bigr((0,T),L^2(\Omega)\bigr)$-norm. Variable time steps are
used to compensate the singularity of the continuous solution near $t=0$.
