Computing integrals involving the matrix exponential

A new algorithm for computing integrals involving the matrix exponential is given. The method employs diagonal Padé approximation with scaling and squaring. Rigorous truncation error bounds are given and incorporated in a Fortran subroutine. The computational aspects of this program are discussed and compared with existing techniques.