A new algorithm for computing integrals involving the matrix exponential is given. The method employs diagonal Pade 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.