figure(1)
clf
N=50;
M=30;
S=1000;
X=zeros(3,N);
A=[ 0 0 .33; .18 0 0 ; 0 .71 .94]
for k=1:M
x=S*[(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*','MarkerSize',20 )
%axis([-10 10 -10 10 -10 10]);
end
[P,D]=eig(A)