function [lambda,x,err] = powers(A,x,eps,N)
%
%  Applies the power method to find the dominant eigenvalue lambda
%  and corresponding eigenvector x of the matrix A.
%  On entry:
%    A         a matrix;
%    x         initial approximation for the eigenvector;
%    eps       accuracy requirement: stops when norm(x-A*x) < eps;
%    N         maximal number of iterations allowed.
%  On return:
%    lambda    modulus of dominant eigenvalue of A;
%    x         eigenvector corresponding with lambda;
%    err       estimate for the errors; 
%              the logarithms of the errors are plotted.
%
%  Example:
%    A = rand(5); x = rand(5,1);
%    [lambda,x,error] = powers(A,x,1.0e-13,50)
%
x0 = x/norm(x);
fprintf('Running the power method...\n');
for i = 1:N
   x = A*x0;
   x = x/norm(x);
   err(i) = norm(x-x0);
   fprintf('  error : %.4e\n', err(i));
   if err(i) < eps
      fprintf('Succeeded in %d steps.\n', i);
      lambda = norm(A*x);
      plot(1:i,log10(err),'*');
      return;
   end;
   x0 = x;
end;
fprintf('Failed to reach accuracy requirement in %d steps.\n', N);
lambda = norm(A*x);
plot(1:N,log10(err),'*');