b=11.677;
Ul=0.5129;
Cea=0.011;
Cel=0.0093;
Cpa=0.0178;
Uamin=0.3;
Uamax=1;
L0=70;  %input('input initial population L0:     ')
P0=30;
A0=65;
n=1000;
jmax=200;
t=zeros(jmax+1,1);
z=zeros(jmax+1,250);
z=zeros(jmax+1,250);
z=zeros(jmax+1,250);
del=(Uamax-Uamin)/jmax;
for j=1:jmax+1
L=zeros(n+1,1);
P=zeros(n+1,1);
A=zeros(n+1,1);
L(1)=L0;
P(1)=P0;
A(1)=A0;
t(j)=(j-1)*del+Uamin;
Ua=t(j);
for i=1:n
L(i+1)=b*A(i)*exp(-Cea*A(i)-Cel*L(i));
P(i+1)=L(i)*(1-Ul);
A(i+1)=P(i)*exp(-Cpa*A(i))+A(i)*(1-Ua);
if (i>750) 
   x(j,i-750)=L(i+1);
   y(j,i-750)=P(i+1);
   z(j,i-750)=A(i+1);
   end
end
end
plot(t,z,'r.','MarkerSize',4)
xlim([Uamin Uamax])
xlabel('\mu_a','FontSize',10), ylabel('Adult population','FontSize',10)
title('Bifurcation diagram for the LPA model')