This paper is concerned with obtaining approximate solutions of fuzzy Fredholm integral equations using Picard iteration method and bivariate Bernstein polynomials. We first present the way to approximate the value of the multiple integral of any fuzzy-valued function based on the two dimensional Bernstein polynomials. Then, it is used to construct the numerical iterative method for finding the approximate solutions of two dimensional fuzzy integral equations. Also, the error analysis and numerical stability of the method are established for such fuzzy integral equations considered here in terms of supplementary Lipschitz condition. Finally, some numerical examples are considered to demonstrate the accuracy and the convergence of the method.