 function [U,T] = schur_proof(A)

% function [U,T] = schur_proof(A)
%
% An algorithmic illustration of our proof of the Schur
% decomposition -- ** do not use in practice! **

 [n n] = size(A);

 if n==1, 
    U = 1;
    T = A;
 else 
    [U,D] = eig(A);
    t11 = D(1,1);
    u1 = U(:,1);
    [Q,R]  = qr(U);    % generate an orthonormal basis for the columns of U
    Q = Q(:,2:end);

    AA = [u1 Q]'*A*[u1 Q];
    t11 = AA(1,1);          % same as D(1,1)
    that = AA(1,2:end)';    % that' = AA(1,2:end)
    Ahat = AA(2:end,2:end);
 
    [Uhat,That] = schur_proof(Ahat);  % recursive call to schur_proof
                                      % with (n-1)-by-(n-1) matrix

    T = [t11 that'*Uhat; zeros(n-1,1) That];
    U = [u1 Q*Uhat];
 end
