Abstract
A numerical solution for neutral delay fractional order partial differential equations involving the Caputo fractional derivative is constructed. In line with this goal, the drift term and the time Caputo fractional derivative are discretized by a finite difference approximation. The energy method is used to investigate the rate of convergence and unconditional stability of the temporal discretization. The interpolation of moving Kriging technique is then used to approximate the space derivative, yielding a meshless numerical formulation. We conclude with some numerical experiments that validate the theoretical findings.