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

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


tic;
xinit = [(1-beta)/2 (1+beta)^2/4];
%xinit = [0.2155    0.3161];
tspan=[0 50];
options=[];
%options=odeset('RelTol',1e-16, 'AbsTol', 1e-16, 'BDF', 'on');
[t,PP] = ode15s('PP',tspan,xinit,options);
toc

figure(3)
hold off
T=25;
N=find(t>T);
plot(t(N:end),PP(N:end,1),'r-',t(N:end),PP(N:end,2),'b-')
%plot(t,x(:,1),'r-',t,x(:,2),'b-')
legend('\itx','\ity','fontsize', 18)
xlabel('Time \it{t}','fontsize', 18)

% Plot solution curve
figure(2)
hold on
T=35;
N=find(t>T);

eps = .3;
[t,PP2] = ode15s('PP',tspan,xinit,options);
N=find(t>T);

plot(PP(N:end,1),PP(N:end,2),'m-',PP2(N:end,1),PP2(N:end,2),'b-')
legend('dd','cc')



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Plot nullclines
%
hold on
%grid off
%caxis([0 1]);
%zmeshsize=.1;
%ameshsize=.5;

% Equilibrium point
%a_star=1+gamma*delta;
%x_star=delta*beta/(a_star*(1-delta));
%z_star=x_star/delta*(1-x_star);
%plot3(x_star,z_star,a_star,'mh','linewidth',5)

%x_star,z_star,a_star;

% x nullcline
x1=0:.05:(1-beta)/2;
y1 =(1-x1).*(beta+x1);
plot(x1,y1,'r--','linewidth',1);
hold on

clear x1 y1
x1=(1-beta)/2:.05:1.1;
y1 =(1-x1).*(beta+x1);
plot(x1,y1,'r-','linewidth',1); 

plot([0 0],[0 beta],'r--',[0 0],[beta 1],'r-','linewidth',2)




% y nullcline
x2=0:.05:1;
y2=(x2./(beta+x2)-delta)./mu;
plot(x2,y2,'b-','linewidth',1);

plot([0 beta*delta/(1-delta)],[0 0],'b-',[beta*delta/(1-delta) 1.1],[0 0],'b--','linewidth',1)



xlabel('x - Prey');ylabel('y - Predator');
%legend('Solution','Prey Nullcline','Predator Nullcline')
axis([0,1.1,0,.7]);


