Abstract: A hybrid numerical method is proposed to evaluate the Green’s functions for displacement of surface or internal arbitrary points in cross-anisotropic multi-layered soil. The wave motion equations in frequency-spatial domain are transformed into state equations in frequency-wavenumber domain using the Fourier transform. Then the dynamic flexibility matrix is obtained. Finally, the inverse Fourier transform is used to evaluate the Green’s functions for displacement of arbitary points in frequency-spatial domain. The proposed algorithm is suitable for arbitrary cross-anisotropic multi-layered soil, and there are no limitations for the number and thickness of the layers. The accuracy of the proposed method is validated by numerical examples, and a parameter analysis for the anisotropy of multi-layered soil is presented.