In this paper, we develop a numerical scheme based on two-dimensional orthonormal Bernstein polynomials (2D-OBPs) to solve two-dimensional nonlinear integral equations of fractional order. The fractional integral considered here is in the Riemann–Liouville sense. By using definition of Riemann–Liouville fractional integral, two-dimensional non- linear fractional integral equations is transformed into two-dimensional nonlinear ordinary integral equations. Operational matrices method based on 2D-OBPs are applied to obtain an approximate solution with high accuracy for these equations. In addition, error analysis of the proposed method is discussed and an upper error bound is provided under weak assumptions. Some linear and nonlinear examples are given to demonstrate the accuracy, efficiency and speed of the suggested method.
