We present an efficient and practical numerical method for solving fractional integro-differential equations with a weakly singular kernel in complex two and three-dimensional domains. The reproducing kernel pseudospectral method is used to discretize the multidimensional regular and irregular spatial domains, while time discretization is done using a third-order backward differentiation formula. We have investigated the unconditional stability of the utilized temporal discretization. Finally, we have used the proposed method to solve the examples presented in the articles and compared the results. The obtained results show the accuracy and ability of the presented method.
