function [l,v,k,error] = pow_it_err(A,tol);
% Usage: [l,v,k] = pow_it(A,tol);
% Computes largest eigenpair (l1,v1) with iteration tolerance tol by
% power iteration.
% Input:
% A = matrix
% tol = tolerance
% Output:
% l = largest eigenvalue
% v = associated eigenvector
% k = number of iterations

[n,n]=size(A);
v = rand(n,1); % Initialize eigenvector
l = (v'*A*v)/(v'*v);
error(1) = norm(A*v-l*v,inf);  % Initialize error
k = 1; maxit = 1000;
while error(k) > tol & k <= maxit
      y = A*v;
      v = y/norm(y,inf);          % Normalize
      l = (v'*A*v)/(v'*v);        % Calculate eigenvalue using Rayleigh quotient 
      k = k+1;
      error(k) = norm(A*v-l*v,inf);
end
if k > maxit
    disp(['Maximum number of iterations reached without achieving tolerance, error = ' sprintf('%0.5g',error)]);
end