A self-consistent nodal method has been developed that directly computes the in-node flux shapes. The method renders the use of an approximation for the transverse leakages no longer necessary. These are obtained directly from the available interface net current shapes, interface flux shapes, and in-node fluxes. The order of the transverse leakage expansion on a set of Legendre polynomials is determined by the order chosen for the method. The results yielded are nearly as accurate (0.02% maximum relative assembly power error) as very fine-mesh benchmark solutions. A comprehensive numerical and analytical analysis of the transverse leakage approximation has been performed. It has been shown that the quadratic leakage approximation can be in error by many times its value. The success of the quadratic leakage approximation is attributed to its small effect on the nodal powers. The theory developed shows that the transverse leakages can have shapes that encompass hyperbolic sines and cosines, and hence that their approximation via quadratic expansions should not always be expected to be adequate. The ILLICO-HO method gives much more information (detailed fluxes and interface currents) than comparable finite difference as well as nodal benchmark solution methods.