%% Parameters

n=200; %Number of time steps. In this case 500 years

%Initializing the starting popluation. We set each initial population to a
%random value at time t=1. This probably makes no sense but we can see what
%happens with the model.


N0(1)=100*rand;
N1(1)=100*rand;
N2(1)=10*rand;
N3(1)=10*rand;
N4(1)=10*rand;
N5(1)=10*rand;
N6(1)=10*rand;
N7(1)=10*rand;
N8(1)=10*rand;
N9(1)=10*rand;
N10(1)=10*rand;
N11(1)=10*rand;
N12(1)=10*rand;


%Setting birth rates. Right now I am assuming most births occur between
%ages 30-35.


b(1)=0;
b(2)=0;
b(3)=.02;
b(4)=.27;
b(5)=1;
b(6)=1;
b(7)=1;
b(8)=.27;
b(9)=.02;
b(10)=.0067;
b(11)=0;
b(12)=0;

%Plotting birth rate data

figure;

naxis=[1,2,3,4,5,6,7,8,9,10,11,12];
plot(naxis,b,'o','LineWidth',2);
set(gca,'FontSize',14)
xlabel('n');
ylabel('b');
axis([1,12,0,2]);
ax = gca;
set(ax,'XTick',[1,2,3,4,5,6,7,8,9,10,11,12]);

%Setting death rates. Right now I am assuming a low mortality rate for
%infants.

d(1)=.5;
d(2)=.1;
d(3)=.1;
d(4)=.05;
d(5)=.05;
d(6)=.05;
d(7)=.1;
d(8)=.1;
d(9)=.2;
d(10)=.2;
d(11)=.3;
d(12)=.3;
d(13)=.4;

figure;
naxis=[0,1,2,3,4,5,6,7,8,9,10,11,12];
plot(naxis,d,'o','LineWidth',2);
set(gca,'FontSize',14)
xlabel('n');
ylabel('d');
axis([0,12,0,1]);
ax = gca;
set(ax,'XTick',[0,1,2,3,4,5,6,7,8,9,10,11,12]);

%% Simulation of the Leslie model

figure;
naxis=[0,1,2,3,4,5,6,7,8,9,10,11,12];
Naxis=[N0(1),N1(1),N2(1),N3(1),N4(1),N5(1),N6(1),N7(1),N8(1),N9(1),N10(1),...
    N11(1),N12(1)];
plot(naxis,Naxis,'o','LineWidth',2);
set(gca,'FontSize',14)
xlabel('n');
ylabel('N');
axis([0,12,0,20]);
ax = gca;
set(ax,'XTick',[0,1,2,3,4,5,6,7,8,9,10,11,12]);

for i=1:n,
    N0(i+1)=N1(i)*b(1)+N2(i)*b(2)+N3(i)*b(3)+N4(i)*b(4)+N5(i)*b(5)+N6(i)*b(6)+...
        N7(i)*b(7)+N8(i)*b(8)+N9(i)*b(9)+N10(i)*b(10)+N11(i)*b(11)+N12(i)*b(12);
    N1(i+1)=(1-d(1))*N0(i);
    N2(i+1)=(1-d(2))*N1(i);
    N3(i+1)=(1-d(3))*N2(i);
    N4(i+1)=(1-d(4))*N3(i);
    N5(i+1)=(1-d(5))*N4(i);
    N6(i+1)=(1-d(6))*N5(i);
    N7(i+1)=(1-d(7))*N6(i);
    N8(i+1)=(1-d(8))*N7(i);
    N9(i+1)=(1-d(9))*N8(i);
    N10(i+1)=(1-d(10))*N9(i);
    N11(i+1)=(1-d(11))*N10(i);
    N12(i+1)=(1-d(12))*N11(i);
    
    Naxis=[N0(i+1),N1(i+1),N2(i+1),N3(i+1),N4(i+1),N5(i+1),N6(i+1),N7(i+1),N8(i+1),N9(i+1),N10(i+1),...
        N11(i+1),N12(i+1)];
    plot(naxis,Naxis,'o','LineWidth',2);
    set(gca,'FontSize',14)
    xlabel('n');
    ylabel('N');
    ax = gca;
    set(ax,'XTick',[0,1,2,3,4,5,6,7,8,9,10,11,12]);
    axis([0,12,0,50]);
    drawnow;
end
