  In this paper we study the finite element approximation of systems of
second-order nonlinear hyperbolic equations. The proposed numerical method
combines a $hp$-version discontinuous Galerkin finite element approximation in
the time direction with an $H^1(\Omega)$-conforming finite element
approximation in the spatial variables. Error bounds at the temporal nodal
points are derived under a weak restriction on the temporal step size in terms
of the spatial mesh size. Numerical experiments are presented to verify the
theoretical results.
