I have a problem where I have $n$ commuting matrices $M_1,\dots,M_n$. It is a well-known fact that commuting matrices are simultaneously diagonalizable/triangularizable. I need to find the eigenvalues of these matrices, but I need to know the eigenvalues grouped up by the common eigenspaces.

In exact arithmetic, this would be as easy as Schur-factorizing $M_1 = UT_1U^*$, and then computing $T_i = U^*M_iU$. However, in floating-point arithmetic, it is my understanding that U may be computed inaccurately when the eigenspaces are poorly conditioned (i.e. when the eigenvalues are clustered).

Is there a stable numerical method for performing this computation?