% Matrix from book---this is a difficult matrix for the QR-eigenvalue algorithm
n = 5;
A = [7.034 -2.271 0 0 0;
    -2.271 2.707 -0.744 0 0;
     0 -0.744 5.804 3.202 0;
     0 0 3.202 -0.464 1.419;
     0 0 0 1.419 -2.082];

%% simplest QR-algorithm 
Ak = A;
for i = 1:50
    [Q,R] = qr(Ak);
    Ak = R*Q
    pause;
end

return;

%% perform QR-algorithm with making matrices smaller
% drop tolerance
tol=1e-3;
% true eigenvalues we want to find
e_real=sort(eig(A));

% maximum iteration nr
kMax = 100;

% initialize
Ak = A;
e_approx = zeros(n,1);
toDo = n;
k = 0;

% loop
while toDo>0
    % find QR-factorization
    [Q,R]=qr(Ak);
    Ak=R*Q;
    Ak
    pause
    % if one eigenvector within tolerance..
    if min(abs(Ak(end,end)-e_real))<tol
        % store eigenvalue and make matrix smaller
        e_approx(toDo)=Ak(end,end);
        toDo = toDo-1;
        Ak = Ak(1:end-1, 1:end-1);
    end
    % if maximal iteration number is obtained, stop
    if k == kMax; break; end
    k = k+1;
end


%% perform QR-algorithm with shift
% Prep
tol=1e-3;
I=eye(n);
% true eigenvalues we want to find
e_real=sort(eig(A));
% maximum iteration nr
kMax=100;
% initialize
Ak=A;
e_approx=zeros(n,1);
toDo=n;
k=0;

% loop
while toDo>0
     % use last element
     % shift=Ak(end,end);
     % Wilkinson
    if size(Ak,1)>=2
        eigBk=eig(Ak(end-1:end, end-1:end));
        [~,IX]=min(abs(eigBk-Ak(end,end)));
        shift=eigBk(IX);
    else
        shift=Ak(end,end);
    end

    [Q,R]=qr(Ak-shift*I);
    Ak=R*Q+shift*I;

    Ak
    pause
    if min(abs(Ak(end,end)-e_real))<tol
        e_approx(toDo)=Ak(end,end);
        toDo=toDo-1;
        Ak=Ak(1:end-1, 1:end-1);
        I=I(1:end-1,1:end-1);
    end

    if k==kMax; break; end
    k=k+1;
end


%% QR collecting Q - i.e. getting the eigenvectors
% Prep
tol=1e-3;

% real eigenvalues we want to find
e_real=sort(eig(A));

% maximum iteration nr
kMax=100;

% initialize
Ak=A;
Qtot=eye(n);

% loop
for k=1:kMax
    [Q,R]=qr(Ak);
    Ak=R*Q;
    Qtot=Qtot*Q;
end

e_approx=diag(Ak);
eVec_approx=Qtot;