Steady-state natural convection heat transfer characteristics from cylinders in a multiply-connected bounded region are clarified. A spectral finite difference scheme (spectral decomposition of the system of partial differential equations, semi-implicit time integration) is applied in numerical analysis, with a boundary-fitted conformal coordinate system through a Jacobian elliptic function with a successive transformation to formulate a system of governing equations in terms of a stream function, vorticity and temperature. Multiplicity of the domain is expressed explicitly.