Using the spherical harmonics method, a universal method for calculating the thermal utilization f is derived by matrix formalism for a heterogeneous reactor with an elementary cell that contains an arbitrary number of concentric cylindrical shells. Matrix elements in the P3 spherical harmonics approximation are obtained and an approximate f-expression by the P3 approximation is derived in a simple and general form as by the ordinary diffusion theory. Finally, f is calculated for two typical cells for a natural uranium heavy-water reactor system and compared with those results for f obtained by the diffusion approximation or the exact P3 approximation.