Optical second harmonic generation (SHG) is a nonlinear optical effect widely used for nonlinear optical microscopy and laser frequency conversion. Closed-form analytical solution of the nonlinear optical responses is essential for evaluating the optical responses of new materials whose optical properties are unknown a priori. A recent open-source code, ♯SHAARP.si, can provide such closed form solutions for crystals with arbitrary symmetries, orientations, and anisotropic properties at a single interface. However, optical components are often in the form of slabs, thin films on substrates, and multilayer heterostructures with multiple reflections of both the fundamental and up to ten different SHG waves at each interface, adding significant complexity. Many approximations have therefore been employed in the existing analytical approaches, such as slowly varying approximation, weak reflection of the nonlinear polarization, transparent medium, high crystallographic symmetry, Kleinman symmetry, easy crystal orientation along a high-symmetry direction, phase matching conditions and negligible interference among nonlinear waves, which may lead to large errors in the reported material properties. To avoid these approximations, we have developed an open-source package named Second Harmonic Analysis of Anisotropic Rotational Polarimetry in Multilayers (♯SHAARP.ml). The reliability and accuracy are established by experimentally benchmarking with both the SHG polarimetry and Maker fringes predicted from the package using standard materials.