This article describes an effective method for the numerical solution of nonlinear time-dependent Volterra-Fredholm integral equations derived from simulating the spatio-temporal spread of an epidemic. At first, the proposed method applies the piecewise linear interpolation technique to discretize the temporal direction. Then, to obtain a full-discrete scheme, the moving least squares (MLS) approach as shape functions in the discrete Galerkin scheme is used to approximate the solution on the domain space. The MLS involves a weighted local least square given on a set of scattered data to estimate multivariate functions. By inheriting the features of the MLS, the offered method can be flexibly employed on non-rectangular domains without any mesh generation on the solution domain. The method’s algorithm is uncomplicated and direct, making it effortlessly executable on a standard PC with regular specifications. The error analysis and convergence rate of the proposed scheme are also discussed. The efficiency and a