A new method for solving the neutron diffusion equations in multiconnected regions with arbitrarily shaped boundaries has been developed by using a compound boundary-fitted coordinate transformation. In the compound boundary-fitted coordinate transformation, inner regions and an outer region in the physical plane are transformed by different coordinate systems. The neutron diffusion equations obtained by the coordinate transformation are solved in the rectangular coordinate system for the outer region, and in the cylindrical coordinate system for the inner regions, so that the boundary conditions are represented accurately and detailed calculations in a particular region can be performed without increasing the number of grid points in other regions.