We develop a finite-difference frequency-domain (FDFD) method to simulate the 2-D SHTE-mode seismoelectric waves. The method approximates the spatial derivatives and non-derivative terms of the wave equations using the weighted-averaging finite-difference operators in a 25-point computational stencil. To suppress the reflections of the seismic and EM waves from the truncated boundaries, we apply the perfectly matched layers (PMLs) surrounding the interior modelled area to absorb the seismic waves, and add several additional layers out of the PMLs to absorb the EM waves. We validate the FDFD method by comparing the FDFD solutions in a two-layer model with the solutions from an analytically-based method. The results show that the FDFD solutions agree excellently with analytical solutions in both the seismic and EM signals, proving that the FDFD method is an efficient and powerful tool in modelling the seismoelectric waves. The FDFD method proposed in this paper requires no quasi-static approximation, and thus can be used to accurately model and interpret the seismoelectric responses in a complex stratum.