  We establish uniform error bounds of a finite difference method for the
Klein-Gordon-Zakharov system (KGZ) with a dimensionless parameter $\varepsilon
\in (0,1]$, which is inversely proportional to the acoustic speed. In the
subsonic limit regime, i.e. $0<\varepsilon \ll 1$, the solution propagates
highly oscillatory waves in time and/or rapid outgoing initial layers in space
due to the singular perturbation in the Zakharov equation and/or the
incompatibility of the initial data. Specifically, the solution propagates
waves with $O(\varepsilon)$-wavelength in time and $O(1)$-wavelength in space
as well as outgoing initial layers in space at speed $O(1/\varepsilon)$. This
high oscillation in time and rapid outgoing waves in space of the solution
cause significant burdens in designing numerical methods and establishing error
estimates for KGZ. By adapting an asymptotic consistent formulation, we propose
a uniformly accurate finite difference method and rigorously establish two
independent error bounds at $O(h^2+\tau^2/\varepsilon)$ and
$O(h^2+\tau+\varepsilon)$ with $h$ mesh size and $\tau$ time step. Thus we
obtain a uniform error bound at $O(h^2+\tau)$ for $0<\varepsilon\le 1$. The
main techniques in the analysis include the energy method, cut-off of the
nonlinearity to bound the numerical solution, the integral approximation of the
oscillatory term, and $\varepsilon$-dependent error bounds between the
solutions of KGZ and its limiting model when $\varepsilon\to0^+$. Finally,
numerical results are reported to confirm our error bounds.
