The transport equation in 1D slab geometry and in the two-angle formulation is analysed in the general case of both anisotropic scattering and a heterogeneous medium. The approach is based on the analytical discrete ordinates technique, which is capable of dealing with highly forward-peaked scattering. The algorithm is extended to deal with multi-layered structures by applying a method derived from the domain decomposition technique, where each single slab is separately solved and the consistency of interface conditions is iteratively sought. To enhance performance, a convergence acceleration of the solution of the linear system is proposed and a study of the dependence of the performance on the adjustable parameters is presented. The optimization of the calculation in terms of computational time and memory requirements is discussed through case studies. The accuracy of the results show great potential as this approach may generate benchmark quality solutions, even for media whose properties change continuously.