A procedure for numerically solving neutron diffusion equations in two-dimensional multiconnected regions with arbitrarily shaped boundaries is developed by using a boundary-fitted curvilinear coordinate transformation. The equation is solved in the transformed rectangular coordinate system where some of the straight coordinate lines represent the boundaries, so that the boundary conditions are represented accurately. Features of the present solution include the automatic generation of coordinate grids, as well as geometric versatility. It can be applied to nuclear design calculations for all types of reactors.