A complete and comprehensive theory of the propagation of second-harmonic (SH) field generated within arbitrary multilayer systems is presented. The theory is based on the radiation of point dipoles whose strength is determined by the distribution of the fundamental field interrogating the multilayer structure. The theory is applied to the study of magnetic-induced second-harmonic generation (MSHG) where the SH field is generated at the interfaces of the multilayer structure as a consequence of the local symmetry breaking, not only due to the structural asymmetry but also the local magnetization. In comparison with the already existing theory based on radiation of an infinitesimally thin polarization sheet, the approach based on the point-dipole radiation presented in this work is more general since it can be applied to the study of systems whose susceptibility tensors, which describe the nonlinear properties of the media, exhibit arbitrary spatial variations. It is shown that the Fourier transform of such variations is closely related to the angular variation of the observable intensity of the far-field SH. Due to the point-dipole nature of the approach used, the theory presented here can be employed in the analysis of MSHG from systems with (buried) magnetic domains with sizes comparable to the wavelength of the second harmonic light, magnetic nanostructures, and others. In addition to the formalism itself, a number of numerical examples is provided and discussed in detail. It is shown that our results are in agreement with those published in the literature for simpler systems. Further capabilities of the theory are demonstrated on a system containing a buried ferromagnetic layer with periodic magnetic domains.