I thank Gosei Furuya and Dan Stanger for their good
suggestions. I was more interested in pedagogy than I
was intrested in building a efficient tool for
matrix exponentiation. For good tools, I recommend
either Gosei's or Dan's code.
These approaches to the matrix exponential require more
machinery than I'm able to build in a sophomore-level DE
course. Instead, I like to use the partial fraction
decomposition (pfd) of the resolvent; without a great
deal of machinery, it's possible to show that the pfd
of the resolvent has the form [see Tosio Kato, page 40]
(M - z I)^-1 = -sum(Pk / (z - zk)) - sum(Dk / (z - zk)^2) + ...,
where z1, ... are the eigenvalues, the Ps are commuting
projections that sum to the identity, and the Ds are nilpotent.
[This is especially easy when the eigenvalues are distinct.
One needn't build the machinery of the Laplace transform
to do this.] For small rigged matrices, finding the pfd of
the resolvent is not that hard to do by hand and a bit of matrix
arithmetic verifies the correctness of the Ps and Ds
(important for hand calculation!)
Although I've done such calculations by hand
for a long time, I'd never tried to build a tool
to do it in general. Things that are simple for
rigged problems, aren't always easy for non-rigged
problems. I played with this for a day, and eventually, I
dumped the Maxima residue code and wrote my own
that works on rational functions only---I was
partially successful; I may or may not try to fix it.
I like the pfd approach to matrix exponentiation for
educational purposes---it's a nice catapult to other
ideas. It works for me.
Of course, the Maxima residue function could use some
work.