function [vec,lastrot,rotangs] = eigpics(A,n)

% This function starts with a random seed vector and repeatedly hits it (n times) with matrix A, and at each step normalizes that vector. If there is a largest-magnitude real eigenvalue the output will approach a unit eigenvector and the rotation angle at each step will approach either 0 or pi, depending on whether the eigenvalue is positive or negative. If the largest magnitude eigenvalue is complex more complicated behavior is to be expected.

% If the matrix is 2x2 or 3x3 this function displays the list of output vectors.

% If n is not specified, its default value is 20.

[a,b]=size(A);

if a~=b ,error('matrix must be 2 by 2 or 3 by 3'),end
if norm(A)==0 , error('matrix must be nonzero'),end

if nargin<2,n=20; end

if a==3

x=rand(3,1);
xold=x/norm(x);
for k=1:n
x=A*x;
x=x/norm(x);
X(k)=x(1);
Y(k)=x(2);
Z(k)=x(3);
rotangs(k)=acos(xold'*x);
xold=x;
end
quiver3(zeros(1,n),zeros(1,n),zeros(1,n),X,Y,Z)
title('Unit Vector Output from Random Seed Vector','FontSize',20,'FontWeight','bold')
vec=x;
lastrot=rotangs(n);

elseif a==2

x=rand(2,1);
xold=x/norm(x);
for k=1:n
x=A*x;
x=x/norm(x);
X(k)=x(1);
Y(k)=x(2);
rotangs(k)=acos(xold'*x);
xold=x;
end
quiver(zeros(1,n),zeros(1,n),X,Y)
title('Unit Vector Output from Random Seed Vector','FontSize',20,'FontWeight','bold')
vec=x;
lastrot=rotangs(n);

else

x=rand(a,1);
xold=x/norm(x);
for k=1:n
x=A*x;
x=x/norm(x);
rotangs(k)=acos(xold'*x);
xold=x;
end
vec=x;
lastrot=rotangs(n);

end
end