A=[-2 1; 1 -2]

fprintf(1,'eigenvalues: %23.15e\n', eig(A));

I2 = eye(2);

v = zeros(2,20);
l = zeros(1,20);
v(:,1) = [1;-5];
l(1) = v(:,1)'*(A*v(:,1))/(v(:,1)'*v(:,1));
fprintf(1,'iteration = %3d, Rayleigh quotient = %23.15e\n', 1, l(1));
pause


for k = 2:5,
  v(:,k) = (A-l(k-1)*I2)\v(:,k-1);
  v(:,k) = v(:,k)/norm(v(:,k));
  l(k) = v(:,k)'*(A*v(:,k));
  fprintf(1,'iteration = %3d, Rayleigh quotient = %23.15e\n', k, l(k));
  pause
end

keyboard
fprintf(1,'Rayleigh quotient = %23.15e\n',l(5));
fprintf(1,'norm(A*v(:,5)-l(5)*v(:,5)) = %23.15e\n', ...
  norm(A*v(:,5)-l(5)*v(:,5))); % residual norm

keyboard

A = [2 1 1;1 3 1;1 1 4]

fprintf(1,'eigenvalues: %23.15e\n', eig(A));

I3 = eye(3);

v = zeros(3,6);
l = zeros(1,6);
v(:,1) = [1;1;1]/sqrt(3);
l(1) = v(:,1)'*(A*v(:,1))/(v(:,1)'*v(:,1));
fprintf(1,'iteration = %3d, Rayleigh quotient = %23.15e\n', 1, l(1));
pause

for k = 2:4,
  v(:,k) = (A-l(k-1)*I3)\v(:,k-1);
  v(:,k) = v(:,k)/norm(v(:,k));
  l(k) = v(:,k)'*(A*v(:,k));
  fprintf(1,'iteration = %3d, Rayleigh quotient = %23.15e\n', k, l(k));
  pause
end

keyboard
fprintf(1,'Rayleigh quotient = %23.15e\n',l(5));
fprintf(1,'norm(A*v(:,5)-l(5)*v(:,5)) = %23.15e\n', ...
  norm(A*v(:,5)-l(5)*v(:,5))); % residual norm