In this article, we introduce a numerical method that combines local radial basis functions partition of unity with backward differentiation formula to efficiently solve linear and nonlinear fractional integro-differential equations on two-dimensional regular and irregular domains. We derive the time-discretized formulation using the backward difference formula. The meshless radial basis function method, particularly the radial basis function partition of unity method, offers advantages such as flexibility, accuracy, ease of implementation, adaptive refinement, and efficient parallelization. We apply the radial basis function partition of unity method to spatially discretize the problem using the scaled Lagrangian form of polyharmonic splines as approximation bases. Numerical simulations demonstrate the efficacy of our method in solving linear and nonlinear fractional integro-differential equations with complex domains and smooth and nonsmooth initial conditions. Comparative analysis confirms the superior performance of our proposed method.
