An example with simple eigenvalues

using JointDiag, LinearAlgebra

N = 4
E = randn(N,N)

Xi = [1 1 1 1   ;
	  1 -1 2 2  ;
	  -1 -1 0 3 ]
	  
M00 = E*diagm(Xi[1,:])*inv(E)
M01 = E*diagm(Xi[2,:])*inv(E)
M02 = E*diagm(Xi[3,:])*inv(E)

M0 = [M00, M01, M02]

X0, E0, Slv = joint_diag(M0, NewtonJointDiag())

The joint eigenvalue vectors are:

julia> X0
3×4 Matrix{Float64}:
 1.0   1.0           1.0   1.0
 2.0   2.0          -1.0   1.0
 3.0  -3.03672e-17  -1.0  -1.0

The common eigenvectors are:

julia> E0
4×4 Matrix{Float64}:
 -0.42709    0.785318     0.18       -0.519269
 -0.824614  -0.45775     -0.0883551  -0.482747
 -0.229514   0.00303585   0.431565    0.588898
 -0.291424  -0.41681      0.879514    0.387962

We verify that the residual error of the decomposition is small:

julia> norm.(M0-[E0*diagm(X0[i,:])*inv(E0) for i in 1:length(M0)])
3-element Vector{Float64}:
 6.926513892003204e-16
 1.0466060440537818e-14
 3.224070268441243e-15