(I’m less worried about using expm1 here and more worried that you are using diagonalization (eigenvectors) to compute matrix exponentials of a non-Hermitian matrix. This can be catastrophically inaccurate if the matrix happens to be nearly defective.)
Note that StaticArrays.jl doesn’t implement a fast eigen the non-Hermitian case, but it does have a fast exp. Maybe compute E = exp(A) and then compute E^n for each n by iterative multiplication?