The fractional reaction–diffusion equation has been used in many real-world applications in fields such as physics, biology, and chemistry. Motivated by the huge application of fractional reaction–diffusion, we propose a numerical scheme to solve the fractional reaction–diffusion equation under different kernels. Our method can be particularly employed for singular and non-singular kernels, such as the Riemann–Liouville, Caputo, Fabrizio–Caputo, and Atangana–Baleanu operators. Moreover, we obtained general inequalities that guarantee that the stability condition depends explicitly on the kernel. As an implementation of the method, we numerically solved the diffusion equation under the power-law and exponential kernels. For the power-law kernel, we solved by considering fractional time, space, and both operators. In another example, we considered the exponential kernel acting on the time derivative and compared the numerical results with the analytical ones. Our results showed that the numerical procedure developed in this work can be employed to solve fractional differential equations considering different kernels.