This paper deals with the numerical solution of nonlinear time-fractional reaction–subdiffusion initial-boundary value problems posed on the space-time domain \(\Omega \times [0,T]\) ; the time derivative is a Caputo fractional derivative of order \(\alpha \in (0,1)\) . First, the problem is transformed into an equivalent integro-differential equation of Volterra type. For this problem a second-order implicit–explicit (IMEX) time-stepping method is combined with the standard 3-point discretisation of the spatial derivative to compute a numerical solution. Under some reasonable assumptions on the data, it is shown that the solution of this method is \(O(N^{-2})\) convergent for \(\alpha \in (1/2,1)\) in the discrete \(L^\infty (0,T; L^2(\Omega ))\) norm when implemented on suitably graded temporal meshes with N points; this result is a significant improvement on the \(O(N^{-(2-\alpha )})\) optimal convergence rate of the well-known L1 scheme. To derive this convergence result a novel error analysis is used, based on subtracting two consecutive steps of our method. Numerical examples are given to illustrate our theoretical result and to provide a comparison with the L1 discretisation. These examples show that \(O(N^{-2})\) accuracy is obtained when \(\alpha \in (0,1]\) ; to prove this result when \(\alpha \in (0,/1/2]\) remains an open problem.
