The radiation exchange in both convex and non-convex enclosures of diffuse gray surfaces is given in the form of a Fredholm boundary integral equation of the second kind. A boundary element method which is based on the Galerkin discretization schem is implemented for this integral equation. Four iterative methods are used to solve the linear system of equations resulted from the Galerkin discretization scheme. A comparison is drawn between these methods. Theoretical error estimates for the Galerkin method has shown to be in a good agreement with numerical experiments.