T = eand it turns out that e^{Mt}a

e^{M} = Σ_{j≥0}M^{j}/j!

e^{Mt} is a square matrix the same size as M.

Some code that agrees with the naïve code.
(Routine `expm` exponentiates a matrix.)

If x is negative the Taylor series for e^{x} is expensive.
If exp(x) directly applies the Taylor series then 1.0/exp(-x) is faster than exp(x) for negative x.
Similarly if you suspect that most of the eigenvalues are negative
inv(expm(−M)) is cheaper than expm(x) to compute.

This hack avoids limitations on M such as symmetry or even lop-sidedness. M can even be complex.