clear all
d=0.7
A=[1.5*(1-d) 0.5*d 0.1*d; 0.75*d 1-d 0.1*d; 0.75*d 0.5*d 0.2*(1-d)]
x=[20; 200; 800]
A*x
y=zeros(3,20);    %build a matrix to store the population data
t=zeros(1,20);    %build a matrix for total population
y(:,1)=x;  %initialize
t(1,1)=sum(y(:,1)); %initialize the sum
for n=1:19
    y(:,n+1)=A*y(:,n);   %update with matrix model
    t(1,n+1)=sum(y(:,n+1));  %update the sume
end
plot([1:20],y(1,:),'*r',[1:20],y(2,:),'xb',[1:20],y(3,:),'+y',[1:20],t(1,:),'og')
legend('paradise','purgatory','inferno','total')

%p=log(t(1,:))
%plot([1:20],p,'*r')

%z=zeros(3,20);
%for m=1:3
%z(m,:)=y(m,:)./t(1,:)
%end
%plot([1:20],z(1,:),'*r',[1:20],z(2,:),'xb',[1:20],z(3,:),'+y')