The Radiative Transfer Equation (RTE) is a multi-dimensional integro - differential equation. It is difficult to obtain an exact analytical solution to the RTE even for simple one - dimensional cases due of its directional nature. As a result, approximate numerical methods must be used to solve the RTE. The two most popular methods that are currently used to solve the RTE are the Method of Spherical Harmonics, PN approximation, and the Discrete Ordinates Method, DOM or SN approximation. However, neither of these methods exhibit good accuracy over the entire range of optical thickness of practical interest. The PN approximation, although it shows good accuracy for optically thick regimes, it is not accurate for optically thin media and in scenarios in which radiation propagation is strongly directional, such as that of a medium bounded by a combination of hot and cold walls. The SN method, on the other hand, does show promising accuracy over a wide range of optical thickness. However, it suffers from ray effects in optically thin media, resulting in locally unphysical solutions. In optically thick media, the strongly coupled directional equations in DOM, renders the method computationally very expensive for obtaining accurate results. Keeping in mind the advantages of each of the afore-mentioned methods and the regimes in which they are accurate, a new robust and computationally efficient hybrid method that has acceptable accuracy over a wide range of optical thickness is proposed, developed, and demonstrated in this thesis. The philosophy of splitting of radiant intensity as used in the Modified Differential Approximation (MDA) is employed here. This philosophy was originally proposed to remove the shortcomings of the P1 approximation in optically thin media. The radiant intensity is split into two components, namely a “ballistic” (wall - emitted) component, and a “diffuse” (medium - emitted) component. Traditionally, the ballistic component is determined using a combination of view-factor based surface-to-surface exchange relations and ray-tracing algorithms, and the diffuse component is determined by invoking the first order spherical harmonics, P1 approximation). Though this method has been shown to be accurate over a wide range of optical thickness, recent studies have shown that this method is prohibitive both from a memory and computational efficiency standpoint for complex three-dimensional geometry with obstructions. The inefficiency stems from the use of the view-factor based approach for determination of the wall-emitted component. In the new hybrid method proposed here, the “wall” component is computed using the Control Angle Discrete Ordinates Method (CADOM), a variant of the Discrete Ordinates Method (DOM), and the “medium” component is evaluated invoking the P1 approximation. This hybrid SN - PN method was validated for both two-dimensional (2D) and three-dimensional (3D) geometries against benchmark Monte Carlo results for gray media in which the optical thickness was varied over a large range. In all cases, the accuracy of the hybrid method was found to be within a few percent of Monte Carlo results, and comparable to the solutions of the RTE obtained directly using CADOM. As a noteworthy advantage, this method was found to be nearly 100 times more efficient than standalone CADOM in optically thick gray media. To demonstrate the method further, 3D non-gray calculations were performed for a gaseous medium with the wide-band model for spectral properties, and the hybrid method was again found to provide substantial computational gains over standard CADOM.