A new iterative method is proposed for solving the three-dimensional neutron diffusion equation. This method reduces the discretization error in the calculation of neutron leakage from a subregion. In addition, when only one fine-mesh point is located in each subregion, this method becomes the same as a fine-mesh finite-difference approximation method. Therefore, it is easy to compare the results of this method with those of a fine-mesh difference approximation. The computer code for this method can be used for calculating both the collapsed neutron flux and fine-mesh difference approximations. The conditions for the convergence of this iterative technique are introduced as a function of the neutron leakage.