% 
% Logistic-Holling Predator-Prey
%
clear all
global eps beta delta mu

eps = 1;
beta = .2;
delta = .2;
mu = 1; 


%%
%% Bifurcation diagram 
%%
xinit = [(1-beta)/2 (1+beta)^2/4];
tspan=[0 40];
epsa=.05;
epsb=1;
K=40;
epsinc=(epsb-epsa)/K;
xextrema=[];
yextrema=[];
for k=0:K
   eps = epsa+k*epsinc;
   if eps < .04
      options=odeset('RelTol',1e-16, 'AbsTol', 1e-16, 'BDF', 'on');
   else 
      options=[];
   end
   [t,PP]=ode15s('PP',tspan,xinit);
   xinit=PP(end,:); % Adaptive initials for different eps: current eps' orbital end is the next eps' orbital initial.
   N=find(t>20); % Throw away transient state and only keep the last t=10 asymptotic state.
   pp=PP(N(1):end,:);
   xextrema=[xextrema;[min(pp(:,1)),max(pp(:,1))]]; % max and min of x
   yextrema=[yextrema;[min(pp(:,2)),max(pp(:,2))]]; % max and min of y
end
epsbif=epsa:epsinc:epsb;
figure(1)
hold off
plot(epsbif',yextrema(:,1),'r-',epsbif',yextrema(:,2),'r-','linewidth',3)
hold on; 
plot([epsb epsb+.3], yextrema(end,:),'r-','linewidth',3) % Extend the equilibrium branch
xlabel('\epsilon','fontsize', 20)
ylabel('Range of \ity','fontsize', 18)



