You might be interested in the Maxima function that
finds the matrix spectral representation (see T. Kato,
"Perturbation Theory for Linear Operators"). The
function 'spectral_rep' returns a list of the
form [eigenvalues, projections, nilpotent].
(%i25) load(linearalgebra)$
(%i26) m : matrix([1,5,1],[0,1,9],[0,0,42])$
(%i27) spectral_rep(m);
(%o27) [[42,1],[matrix([0,0,86/1681],[0,0,9/41],[0,0,1]),matrix([1,0,-86
/1681],[0,1,-9/41],[0,0,0])],matrix([0,5,-45/41],[0,0,0],[0,0,0])]
Reconstruct the matrix from its spectral representation.
(%i28) %[1] . %[2] + %[3];
(%o28) matrix([1,5,1],[0,1,9],[0,0,42])
(%i40) m : matrix([1,5,1],[0,1,9],[0,1,42])$
(%i41) sr : spectral_rep(m)$
(%i42) p : sr[2]$
Show that p[1] is a projection.
(%i43) expand(p[1] . p[1] - p[1]);
(%o43) matrix([0,0,0],[0,0,0],[0,0,0])
(%i44) expand(p[1] . p[2]);
(%o44) matrix([0,0,0],[0,0,0],[0,0,0])
(%i45) np : sr[3]$
Show that the last member of sr is nilpotent.
(%i46) expand(np^^3);
(%o46) matrix([0,0,0],[0,0,0],[0,0,0])
Barton