The flux takes the form of an expansion in spherical harmonics, and the twodimensional spherical harmonics equations to threeorder for anisotropic scattering is drived. The equations is solved using finite element method and compared the solutions with the firstorder diffusion equation and Monte Carlo simulation. The resolut shows that the developed threeorder model with high absorb coefficient is able to significantly improve the diffusion solution in circle geometry, and the radiance distribution close to light source is more accurate.