  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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N = 200;
mu_L = .5;
mu_A = .5;
c_el = .01;
c_ea = .01;
c_pa = .01;
for b = 1.5:.5:16
    param = struct('b',b,'mu_L',mu_L,'mu_A',mu_A,'c_el',c_el,'c_ea',c_ea,'c_pa',c_pa);
    %A0 = log(b*(1-mu_L)/mu_A)/(c_ea+c_pa);
    %L0 = b*A0*exp(-c_ea*A0);
    %P0 = L0*(1-mu_L);
    L0 = log(b*(1-mu_L)/mu_A)/(c_ea+c_pa);  % AH
    A0 = L0*(1-mu_L);                       % AH
    P0 = b*A0*exp(-c_ea*A0);                % AH
    x0 = [L0;P0;A0];
    options = optimset('display','off','Jacobian','on');
    xStar = fsolve(@(x) Fhat(x,param),x0,options);
    
    fprintf('b = %3.1f, xStar = [%g, %g, %g], x0 = [%g, %g, %g]\n',b,xStar(1),xStar(2),xStar(3),L0,P0,A0);
    [y,J] = F(xStar,param);
    e = max(abs(eig(J)));
    if 0
    if e<1
       Display(L0,P0,A0,param,N)
    else
       Display(xStar(1),xStar(2),xStar(3),param,N)
    end
    title(sprintf('b = %5.3f  maxeig = %10.6f',b,e),'fontsize',16)
    xlabel(sprintf('L = %10.6f   P = %10.6f   A = %10.6f',xStar),'fontsize',16)
    legend('L','P','A')
    shg
    pause(1)  
    end
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);

      
  function [y,J] = F(x,param)
% Returns in y and J the value of the function F and its Jacobian at x where
L = x(1); 
P = x(2); 
A = x(3);
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;
% Evaluate F...
y = zeros(3,1);
y(1) = b*A*exp(-c_ea*A)*exp(-c_el*L);
y(2) = L*(1-mu_L);
y(3) = P*exp(-c_pa*A)+A*(1-mu_A);
% Evaluate its Jacobian...
J = zeros(3,3);
J(:,1) = [-c_el*b*A*exp(-c_ea*A)*exp(-c_el*L);(1-mu_L);0];
J(:,2) = [ 0 ; 0 ;exp(-c_pa*A)];
J(:,3) = [b*exp(-c_ea*A)*exp(-c_el*L)-c_ea*b*A*exp(-c_ea*A)*exp(-c_el*L);0;...
          -c_pa*P*exp(-c_pa*A)+(1-mu_A)];
      
  function [y,J] = Fhat(x,param)
% Returns in y and J the value of the function Fhat and its Jacobian at x
[y,J] = F(x,param);
y = y-x;
J = J-eye(3,3);

% up to -2 for incorrect plots
% -1 for incorrect labelling of plots
% -2 for not supplying Jacobian to fsolve
% -1 for parameters not as stipulated





