function d = qreigval(A)
% QREIGVAL  QR algorithm with Wilkinson shifts for the eigenvalues 
%           and deflation.

% Iterate over eigenvalues
for n = length(A):-1:2
  % QR iteration
  s = 0;  I = eye(n);
  while sum( abs(A(n,1:n-1)) ) > eps
    [Q,R] = qr(A-s*eye(n));
    A = R*Q + s*eye(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-A(n,n)));
    s = mu(j);
  end 
  % Deflation
  d(n) = A(n,n);
  A = A(1:n-1,1:n-1);
end
% Last remaining eigenvalue
d(1) = A(1,1);