A time-fractional initial–boundary problem is considered on a bounded spatial domain Ω⊂Rd, where d∈{1,2,3} and Ω is convex or smooth. The differential equation is ∑i=1lqiDtαiu(x,t)−Δu(x,t)=f, where each Dtαi is a Caputo derivative with 0<αl<⋯<α1<1 and the qi are positive constants. A new error analysis is given for the numerical method (L1 scheme in time, finite elements in space) of Huang and Stynes (2020). The error bounds that are derived here remain valid (i.e., are “α-robust”) if α1→1−, unlike the bounds in Huang and Stynes (2020), which blow up as α1→1−. Thus these new error bounds are a more natural and desirable result, since the unknown solution is well behaved as α1→1−.