<p>Fractional diffusion models are widely used to describe anomalous transport phenomena but pose significant computational challenges due to their nonlocal temporal operators. In this work, we propose a high-order iterative grouping algorithm for the efficient numerical simulation of multi-dimensional fractional diffusion equations. The method combines a fourth-order finite difference discretization with a structured grouping strategy that reduces iteration counts and computational overhead while preserving numerical stability. The proposed algorithm is analyzed in terms of stability and convergence and is evaluated through a set of numerical experiments that assess accuracy, runtime, and computational efficiency on refined spatial-temporal grids. Performance results demonstrate that the grouping strategy significantly accelerates convergence compared to standard pointwise iterative methods, particularly for moderately refined grids and long-time simulations. The algorithm demonstrates strong computational efficiency on a single computing device. Due to its block-structured formulation, it has potential for parallel and distributed implementations; however, such HPC-oriented performance is not evaluated in the present study. The nonlocal temporal structure of fractional diffusion equations leads to high computational and memory demands for fine spatial-temporal discretizations, motivating the use of HPC for large-scale simulations.</p>