In this paper, we develop an efficient matrix method based on two-dimensional orthonormal Bernstein polynomials (2D-OBPs) to provide approximate solution of linear and nonlinear weakly singular partial integro-differential equations (PIDEs). First, we approximate all functions involved in the considerable problem via 2D-OBPs. Then, by using the operational matrices of integration, differentiation, and product, the solution of Volterra singular PIDEs is transformed to the solution of a linear or nonlinear system of algebraic equations which can be solved via some suitable numerical methods. With a small number of bases, we can find a reasonable approximate solution. Moreover,we establish some useful theorems for discussing convergence analysis and obtaining an error estimate associated with the proposed method. Finally, we solve some illustrative examples by employing the presented method to show the validity, efficiency, high accuracy, and applicability of the proposed technique.
