A=[-2 1; 1 -2]

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

v = zeros(2,21);
v(1:2,1) = [1;-5];           % starting vector
l(1) = norm(v(:,1),'inf');
v(:,1) = v(:,1)/l(1);       % normalize

for k = 2:20,               % 19 steps of power method
  v(:,k) = A*v(:,k-1);
  l(k) = norm(v(:,k),'inf');
  v(:,k) = v(:,k)/l(k);
  fprintf(1,'iteration = %3d, approx. abs. eig. val. = %23.15e\n', k-1, l(k));
  pause
end

k = 21;                     % step 20
v(:,k) = A*v(:,k-1);
[maxv,idx] = max(abs(v(:,21)));    % find index max abs coefficient
lambda = v(idx,21)/v(idx,20);      % estimate eigenvalue
keyboard
fprintf(1,'approx. eig. val. = %23.15e\n', lambda); % eigenvalue
fprintf(1,'norm(A*v(:,20)-lambda*v(:,20)) = %23.15e\n', ...
  norm(v(:,21)-lambda*v(:,20))); % residual norm

keyboard

A = [2 1 1;1 3 1;1 1 4]
fprintf(1,'eigenvalues: %23.15e\n', eig(A));

v = zeros(3,20);
v(1:3,1) = [1;1;1];         % starting vector
l(1) = norm(v(:,1),'inf');
v(:,1) = v(:,1)/l(1);       % normalize

for k = 2:20,               % 19 steps of power method
  v(:,k) = A*v(:,k-1);
  l(k) = norm(v(:,k),'inf');
  v(:,k) = v(:,k)/l(k);
  fprintf(1,'iteration = %3d, approx. abs. eig. val. = %23.15e\n', k-1, l(k));
  pause
end

k = 21;                     % step 20
v(:,k) = A*v(:,k-1);
[maxv,idx] = max(abs(v(:,21))); % find index max abs coefficient
lambda = v(idx,21)/v(idx,20);   % estimate eigenvalue
keyboard
fprintf(1,'approx. eig. val. = %23.15e\n', lambda); % eigenvalue
fprintf(1,'norm(A*v(:,20)-lambda*v(:,20)) = %23.15e\n', ...
  norm(v(:,21)-lambda*v(:,20))); % residual norm