A=[0 0 .42; .6 0 0; 0 .75 .95] eigs=eig(A) [V,d]=eig(A) 100*V(:,3)/V(3,3) norm(d(1,1)) N = 50; x = [0;0;100]; B = [x]; for n = 1:N-1, x=A*x/eigs(3); % eigs(3) is the real eigenvalue B = [B x]; end title('Buffalo Populations'); xlabel('calf,yearling,adult'); ylabel('relative proportion'); plot(B) clg pause(.1) title('Buffalo orbit'); xlabel('calf'); ylabel('yearling'); plot(B(1,:),B(2,:)) clg pause(.1)