MatrixExp(mat, t = 1, n = 20, k = 3)mat$$E = U exp(D) U^{-1}$$
where $D$ is a diagonal matrix with the eigenvalues of $M$ on the diagonal, $exp(D)$ is a diagonal matrix with the exponentiated eigenvalues of $M$ on the diagonal, and $U$ is a matrix whose columns are the eigenvectors of $M$.
However, if $M$ has repeated eigenvalues, then its eigenvector matrix is non-invertible. In this case, the matrix exponential is calculated using a power series approximation,
$$exp(M) = I + M + M^2/2 + M^3 / 3! + M^4 / 4! + ...$$
For a continuous-time homogeneous Markov process with transition
intensity matrix
$Q$, the probability of occupying state $s$ at time $u + t$
conditional on occupying state $r$ at time $u$ is given by the
$(r,s)$ entry of the matrix $exp(tQ)$.
The series approximation method was adapted from the corresponding
function in Jim Lindsey's R package rmutil.