  We propose a new numerical method to solve linear ordinary differential
equations of the type $\frac{\partial u}{\partial t}(t,\varepsilon) =
A(\varepsilon) \, u(t,\varepsilon)$, where
$A:\mathbb{C}\rightarrow\mathbb{C}^{n\times n}$ is a matrix polynomial with
large and sparse matrix coefficients. The algorithm computes an explicit
parameterization of approximations of $u(t,\varepsilon)$ such that
approximations for many different values of $\varepsilon$ and $t$ can be
obtained with a very small additional computational effort. The derivation of
the algorithm is based on a reformulation of the parameterization as a linear
parameter-free ordinary differential equation and on approximating the product
of the matrix exponential and a vector with a Krylov method. The Krylov
approximation is generated with Arnoldi's method and the structure of the
coefficient matrix turns out to have an independence on the truncation
parameter so that it can also be interpreted as Arnoldi's method applied to an
infinite dimensional matrix. We prove the superlinear convergence of the
algorithm and provide a posteriori error estimates to be used as termination
criteria. The behavior of the algorithm is illustrated with examples stemming
from spatial discretizations of partial differential equations.
