function gamma = qrshift(A)
% QRSHIFT   QR iteration with Wilkinson shift for eigenvalues.
% (demo only--not for practical efficiency)
% Input:
%   A         square matrix
% Output: 
%   gamma     sequence of approximations to one eigenvalue (vector)

n = length(A);
s = 0;  I = eye(n);
k = 1;
while norm(A(n,1:n-1)) > sqrt(n)*eps
  [Q,R] = qr(A-s*I);
  A = R*Q + s*I;
  gamma(k) = A(n,n);
  b = -(A(n-1,n-1)+A(n,n));
  c = A(n-1,n-1)*A(n,n) - A(n-1,n)*A(n,n-1);
  mu = roots([1 b c]);
  [temp,j] = min(abs(mu-gamma(k)));
  s = mu(j);
  k = k+1;
end