This paper introduces a novel numerical scheme for solving the two-dimensional fractional space–time partial differential equation governing anomalous diffusion in the absence of magnetic field gradients. The generalized equation incorporates fractional derivatives in both time (order \(\gamma\) ) and space (order \(\theta\) ), along with a reaction/decay term, providing a comprehensive model for complex diffusion phenomena. For temporal discretization, we employ a high-order time stepping finite difference scheme, ensuring high accuracy in capturing memory effects. For spatial approximation, we propose a new class of shifted Chelyshkov polynomials (in case two dimensional) along to the operational matrices corresponding to left/right Riemann–Liouville and Riesz fractional derivative operators, uniquely constructed for this problem, enabling efficient and stable spectral convergence. The high computational cost of the proposed method, stemming from dense fractional derivative matrices, requires high-performance computing. Its inherent Kronecker structure enables efficient distributed-memory parallelization for large-scale simulations and parameter space exploration. A rigorous convergence analysis is presented for the time semi-discrete scheme, proving that the method achieves unconditional stability and convergence rate of order \(\mathcal{O}(\delta ^{3-\gamma })\) under mild regularity conditions with the time step \(\delta\) . Moreover, numerical experiments demonstrate the superiority of the proposed approach compared to existing techniques, showcasing significantly lower errors and enhanced computational efficiency.
