Theorem (Horn and Johnson, Matrix Algebra); If P is a commuting family of matrices, then there exists a unitary matrix U which simultaneously upper- triangularizes every member of P: .
Observation. The roots of the original system are these associated eigenvalues: , etc. This observation is astonishing, to me; it is essentially a generalization of the companion matrix method for finding roots of univariate polynomials (which is how Matlab finds roots of polynomials).
To see that the eigenvalues give the roots of the polynomial system, let us work in the basis given by the matrix U (for simplicity let's assume all the roots are distinct). Then multiplication of a basis vector by gives us ; multiplication of this by gives ; addition of a constant times gives the value . Now is, in the original basis, the vector of linear coefficients of a polynomial (call it ). We now see that reduces to . If we take , for example one of the Gröbner basis elements, then we know that so it reduces to zero. Therefore , and by the linear independence of the basis we must have that the eigenvalues are roots of f.
(I don't see an easy way to show that all roots are eigenvalues, at this moment).
I don't wish to describe this algebra here, but I refer you to my paper in the proceedings of the ISSAC '95 conference in Montréal for a similar analysis, and to L. González-Vega's paper, and to the references contained therein. We content ourselves here with an example.
(Note that there is no package to do this automatically in Maple. I think that this would make an excellent project, now that I come to think of it).
We end with one last observation, and that is if the roots are distinct, then we can find this matrix U in a reasonably efficient way by finding the eigenvectors of one of these multiplication matrices; the matrix of eigenvectors then gives us our U.
For the example just given, the matrices are
Then the eigenvalues are and , and the desired roots of the system are , , and . The multiple root at zero makes finding the unitary matrix U explicitly quite difficult, but in this case we can sort out which eigenvalue goes with which by simple elimination (because the problem is so small).