clear all;
close all;

T = 50; 
n = 10;
dT = T/n;
a = [ 0.00	0.00024	0.05861	0.28608	0.44791	0.36399	0.22259	0.10457	0.02826	0.00240 ];
b = [ 0.99651	0.99820	0.99802	0.99729	0.99694	0.99621	0.99460	0.99184	0.98700	0 ];

% Leslie matrix
L = [ a ; diag(b(1:n-1)), zeros(n-1,1) ];
% initial age distribution vector
x = ones(n,1);
X(:,1) = x;
t(1) = 1965;
% dynamics of female population
kMax = 20;
for k = 1:kMax
    x = L*x;
    X(:,k+1) = x;
    t(k+1) = 1965 + dT*k;
end
figure(1); plot(t,X); xlabel('years'); ylabel('Number of Canadian women in age classes'); 

% Yield matrix
Y(:,1:kMax-1) = (X(:,2:kMax) - X(:,1:kMax-1))./X(:,1:kMax-1);
figure(2); plot(t(2:kMax),Y); xlabel('years'); ylabel('Yield of Canadian women in age classes'); 

kk = (2040 - 1965)/dT;
s = sum(X(:,kk+1));
z = X(:,kk+1)*100/s;
age = 2.5:5:47.5;
figure(3); plot([age'-2.5,age'+2.5]',[z z]','b'); 
hold on; plot([age'-2.5,age'-2.5]',[zeros(n,1),z]','b');
plot([age'+2.5,age'+2.5]',[zeros(n,1),z]','b'); 
xlabel('age classes'); ylabel('Percentage of the age class in total population');
zNorm = z/z(1); zNorm = zNorm'
