  We develop an algorithm for computing the solution of a large system of
linear ordinary differential equations (ODEs) with polynomial inhomogeneity.
This is equivalent to computing the action of a certain matrix function on the
vector representing the initial condition. The matrix function is a linear
combination of the matrix exponential and other functions related to the
exponential (the so-called phi-functions). Such computations are the major
computational burden in the implementation of exponential integrators, which
can solve general ODEs. Our approach is to compute the action of the matrix
function by constructing a Krylov subspace using Arnoldi or Lanczos iteration
and projecting the function on this subspace. This is combined with
time-stepping to prevent the Krylov subspace from growing too large. The
algorithm is fully adaptive: it varies both the size of the time steps and the
dimension of the Krylov subspace to reach the required accuracy. We implement
this algorithm in the Matlab function phipm and we give instructions on how to
obtain and use this function. Various numerical experiments show that the phipm
function is often significantly more efficient than the state-of-the-art.
