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 L 2 ( 0 , T ) , L 2 ( Ω ) $L^{2}\left ((0,T),L^{2}({\Omega })\right )$ -norm is O ( k 2 − μ 2 + h 2 ) $O(k^{2-\frac {\mu }{2}}+h^{2})$ (that is, short by order μ 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 ∞ ( 0 , T ) , L 2 ( Ω ) $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.