function [mu,x,M,X]=powm91aa(A,x0,Nmax,TOL);
%%function [mu,x,M,X]=powm91aa(A,x0,Nmax,TOL);
%% The power method with Aitken's Delta^2 accerleration
%% 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 (accelerated) approximations to the eigenvalue
%%   X  - matrix of successive approximations to the eigenvectors as rows

k=1;flag=0;x=x0;mu0=0;mu1=0;
I=find(abs(x)==norm(x,inf));
x=x/x(I(1));
while (k<=Nmax & flag==0)
    y=A*x;
    mu =y(I(1));
    if mu-2*mu1+mu0==0;
          flag=1;
    else
          mutil=mu0-((mu0-mu1)^2/(mu-2*mu1+mu0));
    end;
    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;mutil];x=y/y(I(1)); X=[X;x'];mu0=mu1;mu1=mu;
end;
