Abstract. An exponential wave integrator sine pseudospectral method is presented and analyzed for discretizing the Klein-Gordon-Zakharov (KGZ) system with two dimensionless parameters 0 < ε ≤ 1 and 0 < γ ≤ 1 which are inversely proportional to the plasma frequency and the speed of sound, respectively. The main idea in the numerical method is to apply the sine pseudospectral discretization for spatial derivatives followed by using an exponential wave integrator for temporal derivatives in phase space. The method is explicit, symmetric in time, and it is of spectral accuracy in space and second-order accuracy in time for any fixed ε = ε 0 and γ = γ 0 . In the O(1)-plasma frequency and speed of sound regime, i.e., ε = O(1) and γ = O(1), we establish rigorously error estimates for the numerical method in the energy space H 1 × L 2 . We also study numerically the resolution of the method in the simultaneous high-plasma-frequency and subsonic limit regime, i.e., (ε, γ) → 0 under ε γ. In fact, in this singular limit regime, the solution of the KGZ system is highly oscillating in time, i.e., there are propagating waves with wavelength of O(ε 2 ) and O(1) in time and space, respectively. Our extensive numerical results suggest that, in order to compute "correct" solutions in the simultaneous high-plasma-frequency and subsonic limit regime, the meshing strategy (or ε-scalability) is time step τ = O(ε 2 ) and mesh size h = O(1) independent of ε. Finally, we also observe numerically that the method has the property of near conservation of the energy over long time in practical computations.