In this paper, we present an alternative method to Sommerfeld integrals for the calculation of a point source radiation in a multilayered background, based on the method of auxiliary sources, also called fictitious sources method. The method lies upon the decomposition of reflected and transmitted fields on a basis of secondary sources, which amplitudes are determined by applying boundary conditions at each multilayer interface, solving an over-determined system of equations. We present a generalization of the classical FSM to open stratified domains, with arbitrary number of layers and arbitrary position of the point source in the multilayer. Two and three dimensional formalisms are proposed, and compared to exact Sommerfeld formulation. Computation time reduction factors as high as 100 are obtained for 2D problems. Moreover, this method enables the study of non-flat interfaces, which is not possible with Sommerfeld approach. This method can be used for example in order to determine radiation efficiency of emitting devices, or it can be integrated in discrete source type electromagnetic methods for calculation of scattering by a complex target.