During the last several years, a class of algorithms has been developed based on two-dimensional–one-dimensional (2D-1D) decomposition of the reactor transport problem. The current 2D-1D algorithm implemented in the DeCART (Deterministic Core Analysis based on Ray Tracing) code solves a set of coupled 2D planar transport and 1D axial diffusion equations. This method has been successfully applied to several light water reactor analysis problems. However, applications with strong axial heterogeneities have exposed the limitations of the current diffusion solvers used for the axial solution. The work reported in this paper is the implementation of a discrete ordinates (Sn)-based axial solver in DeCART. An Sn solver is chosen to preserve the consistency of the angular discretization between the radial method of characteristics and axial solvers. This paper presents the derivation of the nodal expansion method (NEM)-Sn equations and its implementation in DeCART. The subplane spatial refinement method is introduced to reduce the computational cost and improve the accuracy of the calculations. The NEM-Sn axial solver is tested using the C5G7 benchmark. The DeCART results with the axial diffusion solver shows keff errors of approximately −95, −74, and −110 pcm for the unrodded configuration, rodded configuration A, and rodded configuration B, respectively. These errors decrease to approximately −40, −11, and −12 pcm by using the NEM-Sn solver. In terms of pin power distribution, the use of the NEM-Sn solver has a small effect, except for the heavily rodded configuration. The implementation of the subplane scheme makes it possible to maintain a coarse axial mesh and therefore to reduce the computational cost of the three-dimensional calculations without reducing the accuracy of the solution.