Immediately after defining $e^{\textbf{A}t}$ using the exponential series, chapter 16 of Shaum's Outlines: Differential Equations (3rd Edition) gives a more efficient method for actually computing $e^{\textbf{A}t}$ for any square matrix $\textbf{A}$. However, the book does not give any hint about why the method works, aside from a throwaway remark that it follows "with some effort" from the Cayley-Hamilton theorem that the infinite series for $e^{\textbf{A}t}$ is equivalent to a polynomial in $t$.
If this algorithm had a distinctive name, it would be easy to search the web for more information, but I am not aware if it has a name or not. Here is the algorithm:
Assume that the desired matrix can be written as $e^{\textbf{A}t} = \alpha_{n-1}\textbf{A}^{n-1}t^{n-1} + \alpha_{n-2}\textbf{A}^{n-2}t^{n-2} + \ldots + \alpha_2\textbf{A}^2t^2 + \alpha_1\textbf{A}t + \alpha_0\textbf{I}$, where each $\alpha_i$ is some function of $t$.
Define $r(\lambda) = \alpha_{n-1}\lambda^{n-1} + \ldots + \alpha_2\lambda^2 + \alpha_1\lambda + \alpha_0$. Then, for any eigenvalue $\lambda_i$ of $\textbf{A}t$, $e^{\lambda_i} = r(\lambda_i)$. (Why? I haven't the foggiest idea.)
Find all the eigenvalues of $\textbf{A}t$. Substitute them one by one into the equation in (2) to get $n$ equations in the $\alpha$'s. Solve this system of equations to get the $\alpha$'s. Then substitute them into the equation in (1) and solve for $e^{\textbf{A}t}$.
Step 3 is trivial. However, it's not at all clear why $e^{\textbf{A}t}$ can always be written in form (1), and even if it can, it is even more baffling why the relation in (2) should be true.
Can anyone explain?
It's best to regard $B=At$ as a whole (or assume $t=1$ as shown in the comment), and consider the problem of computing $$e^{B}=\sum_{m=0}^\infty\frac{1}{m!}B^m$$
Because of Cayley-Hamilton, we know that $B^m$ can be written as a polynomial of order $<n$. Hence we may assume $$e^B = \sum_{m=1}^{n-1}\alpha_m(B)B^m$$
If $\lambda$ is an eigenvalue of $B$, then take an eigenvector $v$, we got $Bv=\lambda v$, hence $$e^Bv=\sum_{m=0}^\infty \frac{1}{m!}B^mv=\sum \frac{1}{m!}\lambda^nv=e^{\lambda} v$$
$$e^Bv=\sum\alpha_mB^mv=\sum\alpha_m\lambda^mv=(\sum\alpha_m\lambda^m)v$$
By comparison, we get $e^{\lambda}v=(\sum\alpha_m\lambda^m)v$ and due to $v\not=0$, we must have $$e^{\lambda}=\sum_{m=0}^{n-1}\alpha_m\lambda^m$$
In fact, step (3) is problematic, since $B$ may not have $n$ distinct eigenvalues. (If $B$ does, then to diagonalize $B=P^{-1}\Sigma P, e^B=P^{-1}e^{\Sigma}P$ would be faster.)