function [lambda,V,cnt]=aitkeneig(A,X0,epsilon,max1)
% Power Method with Aitken  Process
%Input    - A is an nxn matrix
%            - X is the 1xn starting vector
%            - epsilon is the tolerance
%            - max1 is the maximum number of iterations
%Output - lambda is the dominant eigenvalue
%            - V is the dominant eigenvector
%            - cnt is the number of iteration

%  NUMERICAL METHODS: Matlab Programs
% (c) 2007 by Xie Liling
%  Complementary Software to accompany the textbook:
%  Information and Computing Science: A Laboratory Course

lambda(1)=0;
cnt=0;
err=1;
k=2;
Y0=A*X0;

%Normalize Y
[m j]=max(abs(Y0));
c(1)=m;
X(1,:)=(1/c(1))*Y0;

t=A*X(1,:)';
Y(1,:)=t';
[m j]=max(abs(Y(1,:)));
c(2)=m;
X(2,:)=(1/c(2))*Y(1,:);

while((cnt<=max1)&(err>epsilon))
   t=A*X(k,:)';
   Y(k,:)=t';
   [m j]=max(abs(Y(k,:)));
   c(k+1)=m;
   X(k+1,:)=(1/c(k+1))*Y(k,:);
   k=k+1;      
   
   % Aitken Process
   pc=c(k-2)-(c(k-1)-c(k-2))^2/(c(k)-2*c(k-1)+c(k-2));
   % pV=X(k-2,:)-(X(k-1,:)-X(k-2,:)).^2./(X(k,:)-2*X(k-1,:)+X(k-2,:));
   lambda(k-1)=pc;
   % V=pV;
   V=X(k,:);
   
   %Update X and lambda and check for convergence
   dc=abs(lambda(k-2)-pc);
   % dv=norm(X(k-1,:)-X(k,:));
   % err=max(dc,dv);
   err=dc;
   
   cnt=cnt+1;   
end
V=(1/lambda(end))*Y(end,:);
V=V';