This paper describes the scattering matrix approach to obtain the solution to electromagnetic field quantities in harmonic multi-layer models. Using this approach, the boundary conditions are solved in such way that the maximum size of any matrix used during the computations is independent of the number of regions defined in the problem. As a result, the method is more memory efficient than classical methods used to solve the boundary conditions. Because electromagnetic sources can be located inside the regions of a configuration, the scattering matrix formulation is developed to incorporate these sources into the solving process. The method is applied to a 3D electromagnetic configuration for verification.