function Y = powerit(A,x0,max,tol)
% power iteration to find the largest eigenpair of A
% x0 initial column vector, length 1
% tol is the tolerance bound
% use infinity norm throughout computations

x = x0; error = 1; i = 0;
while (i <= max) & (error > tol)
    w = A*x;
    x = w./norm(w,inf);
    lambda = x'*A*x/(x'*x);
    error = norm(A*x-lambda*x,inf);
    i = i + 1;
    Y(i,:) = [x;lambda;error].';    % real transpose
end
