  function P2()
% Examines a dynamical system that captures the essentials of the flour
% beetle life cycle.
clc
close all
% Number of time steps in the simulation...
N = 200;
% Model parameters...
b = 11.6772;
mu_L = .5129;
c_el = .0093;
c_ea = .0110;
c_pa = .0178;
L0 = 70;
P0 = 30;
A0 = 70;
% Vary mu_A...
for mu_A = .5:.025:.8
    % Package all the model parameters for convenience...
    param = struct('b',b,'mu_L',mu_L,'mu_A',mu_A,'c_el',c_el,'c_ea',c_ea,'c_pa',c_pa);
    % Watch how L, P, and A vary over the simulation...
    Display(L0,P0,A0,param,N)
    title(sprintf('mu_A = %5.3f',mu_A),'fontsize',16)
    legend('L','P','A')
    shg
    pause(1)   
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  function Display(L0,P0,A0,param,N)
% Plots the L, P, and A values for N steps of the simulation starting
% with initial values L0, P0, and A0.
L = zeros(1,N+1); L(1) = L0;
P = zeros(1,N+1); P(1) = P0;
A = zeros(1,N+1); A(1) = A0;
for k=1:N
   [L(k+1),P(k+1),A(k+1)] = TimeStep(L(k),P(k),A(k),param);
end
plot(1:N+1,L,1:N+1,P,1:N+1,A)
        

  function [L_next,P_next,A_next] = TimeStep(L,P,A,param)
% One "time step" in the simulation.
% Decode the model parameters..
b = param.b; 
mu_L = param.mu_L; 
mu_A = param.mu_A; 
c_el = param.c_el; 
c_ea = param.c_ea; 
c_pa = param.c_pa;
% The step...
L_next = b*A*exp(-c_ea*A)*exp(-c_el*L);
P_next = L*(1-mu_L);
A_next = P*exp(-c_pa*A) + A*(1-mu_A);

 