function [lam,x] = ray_quo(A, x0, nsteps,displ)
%
%  call: [lam,x] = ray_quo(A, x0, nsteps,displ)
%
%  Rayleigh quotient iteration with A, starting with vector x0
%  Returns approximation to eigenvalue and (right) eigenvector.
%
x = x0/norm(x0);      
lam = x'*A*x;
I = eye(length(x0));
for i = 1:nsteps
   y = (A - lam*I)\x;  % approximates a (right) e-vector for e-value near lam
   x = y/norm(y);      % normalize
   lam = x'*A*x;       % Rayleigh quotient approximates eigenvalue
   if displ
      lam
      pause
   end
end
