% This program computes a bifurcation diagram for the holling 
% type II predator-prey model, using K as the bifurcation parameter.

rect = [200 80 700 650]; %fix the window size and position
set(0, 'defaultfigureposition',rect);
global time IC r K s c a d
r=2; s=1; c=0.5; a=2; d=0.1;
option=odeset('AbsTol',1e-11,'RelTol',1e-11);
inc=[0.01:1:2200]';
time=[ inc ] ;
limit=[ 2000:1:2200];
for K=0.1:.05 : 1.8,
IC=[0.5 .2];
[t,U]=ode15s('predprey',time,IC);
u1=U(limit,1);
u2=U(limit,2);
xmin=min(u1);
xmax=max(u1);
ymin=min(u2);
ymax=max(u2);

subplot(211),
plot(K, xmin,'*','MarkerSize',2);
plot(K, xmax,'+','MarkerSize',2);
%axis(ax);
hold on
xlim([.1 1.8]); ylim([0 1.6]);
end
title('Bifurcation diagram for the Holling type II predator-prey model','FontSize',12)
ylabel('x','Rotation',0);

subplot(212),
inc=[0.01:1:2200]';
time=[ inc ] ;
limit=[ 2000:1:2200];
for K=0.1:.05 : 1.8,
IC=[0.5 .2];
[t,U]=ode15s('predprey',time,IC);
u1=U(limit,1);
u2=U(limit,2);
xmin=min(u1);
xmax=max(u1);
ymin=min(u2);
ymax=max(u2);
plot(K, ymin,'*','MarkerSize',2);
plot(K, ymax,'+','MarkerSize',2);
%axis(ax);
hold on
xlim([.1 1.8]); ylim([0 4]);
end
xlabel('K');
ylabel('y','Rotation',0);