  We consider the numerical integration of the matrix Hill's equation.
Parametric resonances can appear and this property is of great interest in many
different physical applications. Usually, the Hill's equations originate from a
Hamiltonian function and the fundamental matrix solution is a symplectic
matrix. This is a very important property to be preserved by the numerical
integrators. In this work we present new sixth-and eighth-order symplectic
exponential integrators that are tailored to the Hill's equation. The methods
are based on an efficient symplectic approximation to the exponential of high
dimensional coupled autonomous harmonic oscillators and yield accurate results
for oscillatory problems at a low computational cost. Several numerical
examples illustrate the performance of the new methods.
