We propose a piecewise-linear, time-stepping discontinuous Galerkin method to solve numerically a time fractional diffusion equation involving Caputo derivative of order μ ź (0, 1) with variable coefficients. For the spatial discretization, we apply the standard continuous Galerkin method of total degree ≤ 1 on each spatial mesh elements. Well-posedness of the fully discrete scheme and error analysis will be shown. For a time interval (0, T) and a spatial domain Ω, our analysis suggest that the error in L2(0,T),L2(Ω)$L^{2}\left ((0,T),L^{2}({\Omega })\right )$-norm is O(k2źμ2+h2)$O(k^{2-\frac {\mu }{2}}+h^{2})$ (that is, short by order μ2$\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(k2 + h2) error bound in the stronger Lź(0,T),L2(Ω)$L^{\infty }\left ((0,T),L^{2}({\Omega })\right )$-norm. Variable time steps are used to compensate the singularity of the continuous solution near t = 0.