We propose a new mixed finite element formulation for solving radiation–conduction heat transfer in optically thick anisotropic media. At this optical regime, the integro-differential equations for radiative transfer can be replaced by the simplified PN approximations using an asymptotic analysis. The conductivity is assumed to be nonlinear depending on the temperature along with anisotropic absorption and scattering depending on both the direction and location variables. The simplified PN approximations are enhanced by considering a diffusion tensor capable of describing anisotropic radiative heat transfer. In the present study, we investigate the performance of the unified and mixed formulations combining cubic P3, quadratic P2, and linear P1 finite elements to approximate the temperature in the simplified P3 model. To demonstrate the performance of the proposed methodology, three-dimensional examples of nonlinear radiation–conduction equations in optically thick anisotropic media are presented. The obtained numerical results demonstrate the accuracy and efficiency of the proposed mixed finite element formulation over the conventional unified finite element formulation to accurately solve the simplified P3 equations in anisotropic media.