There is a reasonable algorithm for the joint diagonalization of a set of matrices, real or complex, self-adjoint or not, because the Jacobi method can be easily extended: we have found a trick by which the Givens angles used in the plane rotations of the Jacobi method can be readily obtained.
We make available the following material regarding joint diagonalization.
A Matlab/Octave function implementing the joint approximate diagonalization of complex matrices (it can also process real matrices, but there are simplifications in this case). As a quick test of what the function does, call it with this little test program
A a function implementing the joint approximate diagonalization of real matrices (offers some computational savings over the more general complex case).
The joint diagonalization procedure is described in this reprint of a SIAM article. The reference is
@article{SC-siam, author = "Jean-Fran\c{c}ois Cardoso and Antoine Souloumiac", journal = "{SIAM} J. Mat. Anal. Appl.", title = "Jacobi angles for simultaneous diagonalization", pages = "161--164", volume = "17", number = "1", month = jan, year = {1996}}
We have also described the first-order perturbation of a joint diagonalizer in this technical report which can be referred to as
@techreport{PertDJ, author = "Jean-Fran\c{c}ois Cardoso", institution = "T\'{e}l\'{e}com {P}aris", title = "Perturbation of joint diagonalizers. Ref\# 94D027", year = "1994"}