A spatial discretization of the Riesz fractional nonlinear reaction–diffusion equation by the fractional centered difference scheme leads to a system of ordinary differential equations, in which the resulting coefficient matrix possesses the symmetric block Toeplitz structure. An exponential Runge–Kutta method is employed to solve such a system of ordinary differential equations. In the practical implementation, the product of a block Toeplitz matrix exponential and a vector is calculated by the shift-invert Lanczos method. Meanwhile, the symmetric positive definiteness of the coefficient matrix guarantees the fast approximation by the shift-invert Lanczos method. Numerical results are given to demonstrate the efficiency of the proposed method.