This question already has an answer here:
So let's say there are two square matrices A and B that commute (AB-BA=0). This implies that there must exist some complete set of vectors that are eigenvectors of both A and B.
Is there any way to tell Mathematica to find the eigenvectors of both A and B? At this point, I can get the eigenvectors of A and then find the eigenvectors of B by taking linear combinations of those for A by hand, but this is obviously not feasible for larger matrices.