This article is the first in a series devoted to the development of efficient and accurate computational tools for the design of beam assemblies for boron neutron capture therapy within the framework of discrete ordinates spectral nodal methods of neutron transport theory. We begin our study with a multi-layer representation of an assembly, and we derive a discrete ordinates matrix operator that replaces without spatial truncation error the entire multi-layer domain in neutron transmission computations. With the matrix operator derived here, we compute without further ado the angular distribution of neutrons leaving the multi-layer assembly, avoiding thus the use of general-purpose discrete ordinates codes founded in the discretization and numerical solution of the neutron transport equation over a number of spatial cells and angular directions throughout the domain. We perform numerical experiments with a four-layer model assembly, and we conclude this article with a discussion and directions for further developments.