function N=Leslie(A,N0,T)
% function N=Leslie(A,N0,T)
%
% A is a parameter vector to be passed on to seqfnc.
%
% This function calculates the sequence up to element T
%  for the difference equation N_(t+1)=seqfnc(N_t) with N_0=N0.
% N and N0 are vectors with a component for each stage.

% define the parameters

p1 = A(1);
p2 = A(2);
f2 = A(3);
f3 = A(4);

% set up the vector of times

t=[0:T];

% compute the sequence values by iteration
% "N=0" is used to erase any values previously stored as "N".

N = 0;
for i=1:3
    N(1,i)=N0(i);
end
for i=1:T
    N(i+1,:)=seqfnc(N(i,:));
end

L = N(:,1);
Y = N(:,2);
A = N(:,3);

% plot the points (t,N) as "x".  Set the axis ranges and labels.

clf

subplot(2,3,1)
plot(t,L,'kx')
xlim([0,T])
xlabel('\it{t}','FontSize',12)
ylabel('\it{L}','FontSize',12,'Rotation',0)

subplot(2,3,2)
plot(t,Y,'bx')
xlim([0,T])
xlabel('\it{t}','FontSize',12)
ylabel('\it{Y}','FontSize',12,'Rotation',0)

subplot(2,3,3)
plot(t,A,'rx')
xlim([0,T])
xlabel('\it{t}','FontSize',12)
ylabel('\it{A}','FontSize',12,'Rotation',0)

subplot(2,3,6)
plot(t,A./(Y+A),'rx')
xlim([0,T])
xlabel('\it{t}','FontSize',12)
ylabel('\it{A/(Y+A)}','FontSize',12)

subplot(2,3,4)
plot(t,L./(Y+A),'kx')
xlim([0,T])
xlabel('\it{t}','FontSize',12)
ylabel('\it{L/(Y+A)}','FontSize',12)

% the code that follows is used to compute the new population levels
% the parameters need to have been defined earlier

    function y=seqfnc(N)
        y(1) = f2*N(2)+f3*N(3);
        y(2) = p1*N(1);
        y(3) = p2*N(2);
    end

end