We have developed new heterogeneous 3D second and fourth order staggered-grid finite-difference schemes for modeling seismic wave propagation in the Laplace-Fourier domain. Recent interest in full waveform inversion in the Laplace-Fourier domain has been the motivation for the development for these types of wave-field simulators. Our approach is based on the principles of an integral equation approximation technique for the velocity-stress formulation in the Cartesian coordinate system. The fourth-order scheme is obtained by the combination of integral identities for the two elementary volumes — “small” and “large” around nodes where the wave variables are defined. The final matrix formulation for the fourth or second-order scheme is performed only for velocity components and the resulting linear system is solved by Krylov iterative methods. We have applied these simulators for the investigation of the wave fields of the SEG/EAGE model in the Laplace-Fourier domain along with other test models. Our eventual goal is to embed it in an inversion scheme for joint seismic-electromagnetic imaging.