<p>Several important physical phenomena in engineering and scientific fields are modeled by non-local fractional models. The development of numerical schemes becomes crucial due to the non-availability of the exact solutions to such models. However, the numerical approximation of these models is challenging and imposes several computational constraints. In this paper, we have devised two highly efficient, energy-stable numerical schemes to solve space fractional reaction-diffusion models. Spatial discretization is performed using a fourth-order matrix transform technique having the advantage of straightforward extension to two and higher spatial dimensions. The time-stepping schemes are developed using an exponential time differencing approach based on a third-order real-pole restricted Padé approximation and a third-order Padé(1,2) approximation. Computationally efficient versions of the schemes are constructed using a splitting technique. Algorithms based on these schemes are constructed, allowing easy coding, and implemented to perform several numerical experiments on problems of practical interest such as the enzyme kinetics equation, Fisher’s equation, and the Allen-Cahn equation. The proposed schemes allow the accurate and efficient simulation of these dynamical models. Solution profiles are plotted to demonstrate the effectiveness of these schemes. Convergence results are computed to validate the accuracy, and central processing unit time is recorded to show the computational efficiency.</p>