Numerical solution of time-fractional reaction-diffusion systems