The four-flux model is a method to solve light radiative transfer problems in planar, possibly multilayer structures. The light fluxes are modeled as two collimated and two diffuse beams propagating forwards and backwards perpendicularly to the layer stack. In the present contribution, we develop a four-flux model relying on a matrix formalism to determine the reflectance and transmittance factors of stacks of components by knowing those of each individual component. This model is also extended to generate the bidirectional scattering distribution function (BSDF) of the stack by considering an incoming collimated flux in any direction, and by taking into account the directionality of the diffuse fluxes exiting from the material at the border components of the stack. The model is applied to opaque Lambertian backgrounds with flat or rough interface, for which analytical expressions of the BSDF are obtained.