This paper presents a formulation of the spherical harmonics method for analyzing radiative heat transfer through participating planar media. The proposed approach expand only the radiative intensity into a finite series of Legendre polynomials, while the scattering phase function is directly averaged over incident and scattered radiations without its expansion in series of Legendre polynomials as classical spherical harmonics method assumes. A matrix algorithm formulation is then implemented, which enable computer of spherical harmonics solution for higher order without difficulty. Radiative heat transfer through anisotropic scattering media under diffuse incidence is examined. The Mie theory for spherical particles is used to account scattering phase function of participating isotropic media. Different scattering phase functions, scattering albedo and optical thicknesses are used in the numerical analysis. Excellent agreements on radiative heat fluxes, hemispherical transmittance and reflectance are obtained between spherical harmonics methods with and without expansion of scattering phase function in series of Legendre polynomials. For high angular discretization, computational time requirements of each of the two spherical harmonics methods are comparable, while the present spherical harmonics formulation is time consuming for low order angular discretization. However, for low and high angular discretizations, spherical harmonics with and without scattering phase function expansion in series of Legendre polynomials computational times are extremely small. The results also indicate that for the same angular discretization order, spherical harmonics methods are more accurate than the discrete ordinates method.