In this paper, we develop a Legendre spectral operational matrix method for the numerical solution of two-dimensional Volterra–Fredholm integro-differential equations subject to mixed boundary conditions. The proposed approach transforms the physical domain onto a reference square and approximates the unknown solution using a tensor-product Legendre polynomial expansion. Exact operational matrices for differentiation and lower-limit integration are constructed, allowing the original integro-differential problem to be reduced systematically to a finite-dimensional algebraic system for the spectral coefficients. The formulation provides a unified treatment of differential, Volterra, and Fredholm operators within a single spectral framework and avoids complicated discretizations of multidimensional integral terms. For a specialized linear form of the problem, rigorous convergence estimates are established in both L2 and L∞ norms under suitable regularity assumptions on the coefficients and kernels. The analysis shows that the dominant convergence behavior is governed by the differential operator, while the integral terms contribute only higher-order consistency effects. Several benchmark examples involving both linear and nonlinear two-dimensional integro-differential equations are presented to demonstrate the performance of the proposed method. Numerical results exhibit rapid spectral-type error decay as the polynomial degree increases, with the numerical errors approaching machine precision for moderate truncation orders. These results confirm the accuracy, efficiency, and reliability of the proposed Legendre spectral operational matrix framework for solving a broad class of multidimensional integro-differential equations with nonlocal operators.
Ishtiaq Ali (Tue,) studied this question.