<p>This work deals with a numerical approach based on a nodal Galerkin spectral method, which is a proper combination of the Legendre spectral method and the L1/L1-2 discretizations to address a class of weakly singular integro-partial differential equations in one and two dimensions having non-local temporal operators. Such equations are fundamental in modeling processes that exhibit memory and hereditary properties, often encountered in fields like anomalous diffusion and viscoelasticity. The primary goal of this work is to develop a stable high-order numerical scheme that can efficiently solve these problems. The numerical approach follows by combining the temporal semi-discretization and the spectral method that involves both trial and test functions being chosen as the Lagrange interpolating polynomials defined at the Legendre–Gauss–Lobatto points. In addition, we investigate their effect on numerical estimates. The L1 or L1-2 scheme is used to discretize the time-fractional operator of the semi-discrete problem, whereas to deal with the weakly singular Volterra operator, a modified quadrature rule (which is a combination of the composite trapezoidal approximation and the midpoint rule) is employed. In contrast, the Legendre spectral method is utilized for an efficient approximation in the space directions for 1D and 2D problems, for which the rate of convergence in space increases with respect to the smoothness of the solution. This spectral approximation involves the order of Legendre polynomials. In addition, it is numerically observed that the fully discretized approximation leads to a stable and convergent numerical solution for problems having an initial layer in time, which remains reliable regardless of the chosen time step size. Moreover, this hybrid numerical approach leads to higher-order accuracy in space and time, where the accuracy rate depends on the degree of the polynomial approximation, the regularity of the solution space, and the order of the weakly singular kernel. Theoretical findings are supported by extensive numerical experiments of various complex fractional integro-differential equations having weakly singular kernels with different regularity conditions on the exact solutions, and it shows the strong effectiveness of the present approach.</p>