The even parity formulation (EPF) of the discrete ordinates method (DOM) is used to simulate radiative heat transfer in two-dimensional enclosures containing an absorbing-emitting and scattering medium. The discrete ordinates equations for the EPF are second-order differential equations and they are spatially discretized using a second-order central difference scheme. At the boundary, a higher-order upwind scheme is employed to prevent solution instability and minimize errors. The matrix solver of the discretized equations is based on a preconditioned conjugate gradients method. To investigate the accuracy and efficiency of the EPF of the DOM, several two-dimensional benchmark problems with an absorbing-emitting and scattering medium enclosed by gray walls are considered. By taking an appropriate numerical treatment, the numerical results from the EPF appear to compare favorably with other available solutions. However, the even parity solution usually requires more CPU time and iterations to converge in comparison with the conventional DOM, especially for the case with a small optical thickness.