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) ];

% finding dominant eigenvalue
format long;
[P,D] = eig(L); 
eigenvalues = abs(diag(D))
[lambdaPlus, indPlus] = max(max(D))
zPlus = P(:,indPlus)/P(1,indPlus); 
zPlus = zPlus'

% the iterative power method for dominant eigenvalue
x = ones(n,1);
kMax = 20;
for k = 1:kMax
    x = L*x;
    c = max(x);
    x = x/c;
    fprintf('k = %2.0f, c = %10.8f\n',k,c);
end
Eigenvector = x'