% choose a random (nonnegative) matrix  and starting vector
%% nonnegative guarantees a dominant real eigenvalue
clf
clear all
n=3;
A=10*rand(n);
x=randn(n,1);
[y,i]=max(abs(x));
xold=x/x(i);
xs=[xold]
   if n==2,
	plot(xs,'or')
   elseif n==3,
	plot(xs,'or')
   end
   pause
xnew=A*xold
[y,i]=max(abs(xnew));
xnew=xnew/xnew(i);
xs=[xs xnew]
   if n==2,
	plot(xs,'or')
   elseif n==3,
	plot(xs,'or')
   end
   pause
while norm(xnew-xold) > 1e-8,
   xold=xnew;
   xnew=A*xold;
   [y,i]=max(abs(xnew));
   xnew=xnew/xnew(i)
   xs=[xs xnew]
   if n==2,
	plot(xs,'or')
   elseif n==3,
	plot(xs,'or')
   end
   pause
end
temp=A*xnew;
disp(['eigenvector estimate is ',num2str([xnew'])])
disp(['eigenvalue estimate is ',num2str(temp(i))])
disp('now verify eigenvalue/eigenvector approximation')
A*xnew - temp(i)*xnew

