global mu_l mu_a c_ea c_el c_pa b mydatafit myDays

% Beetles, Cannibalism, and Chaos
% Solution to Problem 2
% Dianne O'Leary 12/2006
% Solve a nonlinear equation to determine the equilibrium population
% for the LPA model with these choices of parameters:

c_el = .01;
c_ea = 0.01;
c_pa = 0.01;
mu_l = .5;
mu_a = .5;

bvals = 1:.5:20;

% The following loop solves for the population at each value of b.

j=0;
for b = bvals,
   j=j+1;

% The initial point Lstar;Pstar;Astar is the equilibrium population
% when c_el = 0.

   c_el = 0;
   Astar = log(b*(1-mu_l)/mu_a)/(c_ea+c_pa);
   Lstar = b*Astar*exp(-c_ea*Astar);
   Pstar = Lstar*(1-mu_l);
   evallpa([Lstar;Pstar;Astar])
   c_el = 0.01;
   x0 = [Lstar;Pstar;Astar];
   xl(:,j) = x0;
   xf(:,j) = fsolve(@evallpa,x0,optimset('Jacobian','on'));
 
% The equilibrium population is stable if all eigenvalues of the
% Jacobian matrix for the function 
%
%   [ Lnew ]      [ L ]
%   [ Pnew ]  = F [ P ]
%   [ Anew ]      [ A ]
%
% lie inside the unit circle.

   [v,Jacob] = evallpa(xf(:,j));
   Jacobeiga = abs(eig(Jacob+eye(3)));

   if (max(Jacobeiga)<1)
      stab(j) = 1;
   else
      stab(j) = 0;
   end

end

% Plot the results: population as a function of b.
% Mark the stable ones with plusses.

figure
plot(bvals,xf(1,:),'b',bvals,xf(2,:),'g',bvals,xf(3,:),'r')
hold on
ind = find(stab);
plot(bvals(ind),xf(1,ind),'b+',bvals(ind),xf(2,ind),'g+', ...
     bvals(ind),xf(3,ind),'r+')
legend('Larvae','Pupae','Adults')
xlabel('b')
ylabel('Population')
title('Equilibrium population as a function of b')