Page 155
L=[0 0 0 0 .08 .28 .42;
   .657 0 0 0 0 0 0;
   0 .930 0 0 0 0 0;
   0 0 .930 0 0 0 0;
   0 0 0 .930 0 0 0;
   0 0 0 0 .935 0 0;
   0 0 0 0 0 .935 0]
L^(10)

[evect,eval]=eig(L)
lambda=eval(1)
pf=evect(:,1)
% get pf=0.8586 and eigenvector=
% [-0.3930 -0.3007 ... -0.4532], normalize
% (mult. by a const.) so leading term is 1
pf=pf/pf(1)

                Page 156
V=[1;
L(2,1)/lambda;
L(2,1)*L(3,2)/lambda^2;
L(2,1)*L(3,2)*L(4,3)/lambda^3; 
L(2,1)*L(3,2)*L(4,3)*L(5,4)/lambda^4;
L(2,1)*L(3,2)*L(4,3)*L(5,4)*L(6,5)/lambda^5;
L(2,1)*L(3,2)*L(4,3)*L(5,4)*L(6,5)*L(7,6)/lambda^6]

p=ones(7,1) % column vector of 1's
(L^10)*p
lambda^10*V


                 Exercises/Experiments
1.
n=0:1:100;
P0=(100-n).*(25+n);
plot(n,P0);
P=P0'; % rows=age, columns=time
% no base 0 indexing so P(n,t)=
%   number aged n-1 in year t-1
mu=.0524*(exp(.03.*n)-1);
%   mu(n) applies to age n-1
plot(n,mu);
for t=1:3
   total(t)=sum(P(:,t));
   P(1,t+1)=(1.9/20)*sum(P(22:31,t));
   for n=2:101
     P(n,t+1)=(1-mu(n-1))*P(n-1,t);
   end; end;
total(1) % starting year
total(4)=sum(P(:,4)) % 3 years later