In the present paper, an attempt was made to develop a numerical method for solving a general form of two-dimensional nonlinear fractional integro-differential equations using operational matrices. Our approach is based on the hybrid of two-dimensional block-pulse functions and two-variable shifted Legendre polynomials. Error bound and convergence analysis of the proposed method are discussed. We prove that the order of convergence of our method is O(1 2 2M−1 N M M!). The presented method is tested by seven test problems to demonstrate the accuracy and computational efficiency of the proposed method and to compare our results with other well-known methods. The comparison highlighted that the proposed method exhibits superior performance than the existing methods, even using a few numbers of bases.