clf
N=20;
M=50;
S=100;
X=zeros(3,N);
A=[.9 .09 .01; .01 .9 .09; .09 .01 .9]
for k=1:M
x=[(rand);(rand);(rand)];
x=S*x/sum(x)
 for j = 1:N
    X(:,j)=A^(j-1)*x;     
 end   
plot3(X(1,:),X(2,:),X(3,:),'g-*')
hold on
plot3(X(1,1),X(2,1),X(3,1),'b*')  
plot3(X(1,N),X(2,N),X(3,N),'r*')
%axis([-10 10 -10 10 -10 10]);
end
[P,D]=eig(A)