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
and
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).