In this paper, we present a new highly efficient numerical algorithm for nonlinear variable-order space fractional reaction–diffusion equations. The algorithm is based on a new method developed by using the Gaussian quadrature pole rational approximation. A splitting technique is used to address the issues related to computational efficiency and the stability of the method. Two linear systems need to be solved using the same real-valued discretization matrix. The stability and convergence of the method are discussed analytically and demonstrated through numerical experiments by solving test problems from the literature. The variable-order diffusion effects on the solution profiles are illustrated through graphs. Finally, numerical experiments demonstrate the superiority of the presented method in terms of computational efficiency, accuracy, and reliability.