The spectral numerical mode-matching (SNMM) method is developed to simulate the 3D layered multi-region structures. The SNMM method is a semi-analytical solver having the properties of dimensionality reduction to reduce computational costs; it is especially useful for microwave and optical integrated circuits where fabrication is often done in a layered structure. Furthermore, at some layer interfaces, very thin surfaces such as good conductor surfaces and metasurfaces can be deposited to achieve desired properties such as high absorbance and/or anomalous reflection/refraction. In this work, the 3D SNMM method is further extended from a single interface to multiple layers so that the electromagnetic propagation and scattering in the longitudinal direction is treated analytically through reflection and transmission matrices by using the eigenmode expansions in the transverse directions. Therefore, it effectively reduces the original 3D problem into a series of 2D eigenvalue problems for periodic structures. We apply this method to characterize metasurfaces and lithography models, and show that the SNMM method is especially efficient when the longitudinal layer thicknesses are large compared with wavelength. Numerical experiments indicate that the SNMM method is highly efficient and accurate for the metasurfaces and the lithography models.Index Terms-Bloch (Floquet) periodic eigenmodes, lithography, metasurface, mixed spectral element method, spectral numerical mode-matching method.