% ============================================
% plot the phase plane for   
%         x'= x(a-bx-cy)
%         y'= y(q-rx-sy)
%  Lotka-Volterramodel of competition
%   see Strogatz section 6.4 (pp.155-159).
% ============================================

global a b c q r s

% set parameters
a=3; b=1; c=2; 
q=2; r=1; s=1;


% compute and plot vector feild
% -----------------------------

% define grid and calculate v-fld on grid
xmin=-0.2; xmax=3.1; dxv=0.2;
ymin=-0.2; ymax=2.1; dyv=0.2;
[x,y] = meshgrid(xmin:dxv:xmax, ymin:dyv:ymax);
dxdt = x.*(a - b*x - c*y);                                      
dydt = y.*(q - r*x - s*y);

% plot v-fld
figure(1)
clf
quiver(x,y,dxdt,dydt,'m');   
hold on
axis([xmin,xmax,ymin,ymax]);
%axis equal
xlabel('x');
ylabel('y');

pause;


% plot nullclines
% ---------------

% x-nullcines
plot([0,0],[ymin,ymax],'k');
xg=[xmin:0.01:xmax];
xn=(a-b*xg)/c;
plot(xg,xn,'k');

% y-nullcines
plot([xmin,xmax],[0,0],'k--');
yn=(q-r*xg)/s;
plot(xg,yn,'k--');
axis([xmin,xmax,ymin,ymax]);

pause;


% solve system od ODEs and plot trajectories
% ------------------------------------------

tspan=[0 10];   %time = 0 to 10
x10=[0.2; 0.1]; % initial conditions (x(0),y(0))
[t1,x1]=ode45(@lvcomp,tspan,x10); % solve system of ODEs in lvcomp.m
plot(x1(:,1),x1(:,2));

pause

% corresponding plot for x and y vs time
figure(2)
clf
plot(t1,x1(:,1))
hold on
plot(t1,x1(:,2),'--')

pause;

% more trajctories
figure(1)

tspan=[0 10];   %time = 0 to 10
x20=[2.0; 1.5]; % initial conditions (x(0),y(0))
[t2,x2]=ode45(@lvcomp,tspan,x20); % solve system of ODEs in lvcomp.m
plot(x2(:,1),x2(:,2));

pause

for x0=0.0:0.1:xmax
    x30=[x0; 0.05]; 
    [t3,x3]=ode45(@lvcomp,tspan,x30); 
    plot(x3(:,1),x3(:,2));
end;

for x0=0.0:0.2:xmax
    x30=[x0; ymax];
    [t3,x3]=ode45(@lvcomp,tspan,x30); 
    plot(x3(:,1),x3(:,2));
end;

for y0=0.0:0.1:ymax
    x30=[0.05; y0]; % initial conditions (x(0),y(0))
    [t3,x3]=ode45(@lvcomp,tspan,x30); 
    plot(x3(:,1),x3(:,2));
end;

for y0=0.0:0.2:ymax
    x30=[xmax; y0]; 
    [t3,x3]=ode45(@lvcomp,tspan,x30); 
    plot(x3(:,1),x3(:,2));
end;




