A fast and accurate approach, based on the fast all-modes method (FAM) and the numerical modified steepest descent path method (NMSP), was previously used to calculate the spatial Green's function for a single-layer lossless dielectric medium over a perfect electric conductor. This paper successfully extends that approach to two new cases. The first is the multilayer case where the medium has an arbitrary number of layers. The second is lossy media over an imperfect conductor. The FAM locates all modes accurately on the complex plane. The modes include surface wave modes, leaky wave modes, and improper modes. For a typical six-layer case over a ground plane, the FAM requires only 2.265 s of pre-processing that includes computing 200 mode locations by using a P4 3.2-GHz PC with Matlab. The NMSP is then used to evaluate the steepest descent path integral. Accuracy within 0.2% is achieved in comparison with the benchmark calculations. Within this context of accuracy, the total CPU per distance point is less than 7.6 ms for distances larger than 0.02 free-space wavelength and is less than 2.7 ms for distances larger than two free-space wavelengths. This method is shown to be fast and accurate, even for large-distance interactions in the multilayer and lossy media.