% construct an orthogonal matrix
>> [q,r]=qr([1, -1, 1, 0; 2, 3, 0, 1; 4, 3, 0, 0; -1, 2, 0, 0])
q =
   -0.2132    0.4706    0.8562         0
   -0.4264   -0.4579    0.1455   -0.7664
   -0.8528   -0.0763   -0.1704    0.4877
    0.2132   -0.7504    0.4655    0.4180
r =
   -4.6904   -3.1980   -0.2132   -0.4264
         0   -3.5739    0.4706   -0.4579
         0         0    0.8562    0.1455
         0         0         0   -0.7664

% create simulataneous diagonalizable matrices a,b
% that are not already diagonal by rotating 
% two diagonal matrices into frame with unit vectors q
>> a=q*diag([5,5,-7,9])*q'
a =
   -3.7972   -1.4947    1.7509   -4.7829
   -1.4947    7.0956   -1.1977   -2.0942
    1.7509   -1.1977    5.6030    1.7675
   -4.7829   -2.0942    1.7675    3.0986
>> b=q*diag([-4,3,3,1])*q'
b =
    2.6818   -0.6364   -1.2727    0.3182
   -0.6364    0.5525   -1.7979    1.2771
   -1.2727   -1.7979   -2.5666    0.8650
    0.3182    1.2771    0.8650    2.3323

% check that the commutator of a,b is zero (they commute)
% note diagonal matrices always commute
>> a*b-b*a
ans =
   1.0e-14 *
         0    0.3553   -0.2665    0.1776
   -0.1776    0.1776         0    0.6217
    0.1776         0         0    0.1776
         0   -0.4885   -0.1776   -0.0888

% diagonalize a using matlab
%  va is a matrix with eigenvectors as rows
%  da is a matrix with eigenvalues along the diagonal
>> [va,da]=eig(a)
va =
    0.8562   -0.0000    0.5166    0.0984
    0.1455    0.7664   -0.2411    0.5208
   -0.1704   -0.4877    0.2824    0.8473
    0.4655   -0.4180   -0.7715   -0.0336
da =
   -7.0000         0         0         0
         0    9.0000         0         0
         0         0    5.0000         0
         0         0         0    5.0000

% check that each eigenvector is multiplied (stretched)
%  by the corresponding eigenvalue
% the ./ is used to divide matrices element-by-element
>> a*va
ans =
   -5.9935   -0.0000    2.5831    0.4922
   -1.0183    6.8977   -1.2055    2.6042
    1.1929   -4.3894    1.4121    4.2365
   -3.2586   -3.7624   -3.8575   -0.1682
>> a*va ./ va
ans =
   -7.0000    3.8383    5.0000    5.0000
   -7.0000    9.0000    5.0000    5.0000
   -7.0000    9.0000    5.0000    5.0000
   -7.0000    9.0000    5.0000    5.0000

% Is va orthogonal?  Almost, except for the two
% degenerate eigenvalues.  Any linear combination of
% them is also an eigenvector.  We could use the
% Gram-Schmidt method to create an orthonormal basis,
% but matlab did bother to do so.
>> va'*va
ans =
    1.0000   -0.0000    0.0000    0.0000
   -0.0000    1.0000    0.0000   -0.0000
    0.0000    0.0000    1.0000    0.1905
    0.0000   -0.0000    0.1905    1.0000

% Verify one more time that  v^(-1) a v  is diagonal
>> va\a*va
ans =
   -7.0000    0.0000    0.0000    0.0000
         0    9.0000    0.0000    0.0000
    0.0000         0    5.0000   -0.0000
    0.0000    0.0000   -0.0000    5.0000

% Can b be diagonalized by the same transformation va?
% No, but mostly.  The last two eigenvectors don't work.
>> bk=va\b*va
bk =
    3.0000   -0.0000   -0.0000   -0.0000
         0    1.0000   -0.0000    0.0000
    0.0000   -0.0000    2.3186   -1.6063
   -0.0000    0.0000   -2.6804   -3.3186

% Another way to see that the first two columns
% are eigenvectors of b as well as a
>> b*va ./ va
ans =
    3.0000    1.3194    1.8079  -11.7486
    3.0000    1.0000    8.1091   -2.5751
    3.0000    1.0000   -5.7228   -3.8540
    3.0000    1.0000    2.2017  -40.1478

% Diagonalize the lower right 2x2 block of bk
>> [vk,dk]=eig(bk(3:4,3:4))
vk =
    0.9206    0.2464
   -0.3905    0.9692
dk =
    3.0000         0
         0   -4.0000

% This is the transformation that diagonalizes all of bk
>> uk=blkdiag(1,1,vk)
uk =
    1.0000         0         0         0
         0    1.0000         0         0
         0         0    0.9206    0.2464
         0         0   -0.3905    0.9692

% Just to show that it does indeed diagonalize bk
>> uk\bk*uk
ans =
    3.0000   -0.0000   -0.0000   -0.0000
         0    1.0000   -0.0000    0.0000
    0.0000   -0.0000    3.0000    0.0000
   -0.0000    0.0000   -0.0000   -4.0000

% Combine the two steps need to diagonalize b
>> w=va*uk
w =
    0.8562   -0.0000    0.4372    0.2227
    0.1455    0.7664   -0.4253    0.4454
   -0.1704   -0.4877   -0.0709    0.8908
    0.4655   -0.4180   -0.6971   -0.2227

% Note that unlike va, the eigenvectors in w
% are orthogonal.  This is because w diagonalizes
% both a and b, and is nondegenerate for the combination.
>> w'*w
ans =
    1.0000   -0.0000         0    0.0000
   -0.0000    1.0000    0.0000   -0.0000
         0    0.0000    0.8630    0.0000
    0.0000   -0.0000    0.0000    1.0910

% Just to check that w diagonalizes b ...
>> w\b*w
ans =
    3.0000   -0.0000    0.0000   -0.0000
    0.0000    1.0000   -0.0000    0.0000
    0.0000   -0.0000    3.0000   -0.0000
   -0.0000    0.0000   -0.0000   -4.0000
% ... and also a !!
>> w\a*w
ans =
   -7.0000    0.0000    0.0000    0.0000
         0    9.0000    0.0000    0.0000
    0.0000   -0.0000    5.0000   -0.0000
    0.0000    0.0000   -0.0000    5.0000
    
% Note the overlap between the degenerate
% eigenspaces of a and b.
