function [mu,x,M,X]=powmet91(A,x0,Nmax,TOL);
%%function [mu,x,M,X]=powmet91(A,x0,Nmax,TOL);
%% The power method for finding the dominant eigenvalue and
%% corresponding eigenvector of the nxn matrix A.
%% INPUTS
%%   A - the given matrix
%%   x0 - a starting vector as a column vector
%%   Nmax - the maximum number of iterations
%%   TOL - tolerence in the stopping criterion
%% OUTPUTS
%%   mu - the approximate eigenvalue
%%   x - the approximate coresponding eigenvector
%%   M  - matrix of successive approximations to the eigenvalue
%%   X  - matrix of successive approximations to the eigenvectors as rows

k=1;flag=0;x=x0;
I=find(abs(x)==norm(x,inf));
x=x/x(I(1));
while (k<=Nmax & flag==0)
    y=A*x;
    mu =y(I(1));
    I=find(abs(y)==norm(y,inf));
    if y(I(1))==0;
       flag=1;
       message='A has eigenvalue 0 with eigenvector'
       x
       message='Select new x0 and restart'
    end;
    err= norm(x-y/y(I(1)),inf);
    if err<TOL
       flag=1;
       message='Tolerence achieved'
    end;
    k=k+1;M=[M;mu];x=y/y(I(1)); X=[X;x'];
end;
