function fig232
clc
clf
% Fig. 2.32
% competing species model
disp('Page 167 Fig2.32: Competing Species Model')
disp('dx/dt = 2*x*(1-x/2)-x*y')
disp('dy/dt = 2*x*(1-x/2)-x.*y')
disp('To be shown: vector fields & equilibrium points')
pause;
%
rmin = 0;
rmax = 3;
fmin = 0;
fmax = 3;
dr = 0.2;
df = 0.2;
[x,y] = meshgrid(rmin:dr:rmax,fmin:df:fmax);
u = 2.*x.*(1.-x/2.)-x.*y;
v = 3.*y.*(1.-y/3.)-2.*x.*y;
leng = sqrt(u.^2+v.^2);
%
% normalize version
quiver(x,y,u./leng,v./leng,0.4)
axis equal tight
xlabel 'x'
ylabel 'y'
title 'Vector Field for competing species model'
%
% solve for the equilibrium points
[x0,y0] = solve('2*x0*(1-x0/2)-x0*y0','3*y0*(1-y0/3)-2*x0*y0=0', 'x0','y0');
[x0 y0]
hold on;
plot(double(x0),double(y0),'go')
%
disp('To be shown: solution curves with initial condition lies in [0,1]')
pause;
% solve for the solution curve
fode = inline('[2*y(1)*(1-y(1)/2)-y(1)*y(2);3*y(2)*(1-y(2)/3)-2*y(1)*y(2)]', 't', 'y');
t0 = 0;
tmax = 10;
subplot(1,2,1);
   quiver(x,y,u./leng,v./leng,0.4)
hold on
plot(double(x0),double(y0),'go')
for a=0.1:0.3:1
    for b = 0:0.4:1
    [t,ysol] = ode45(fode,[t0 tmax], [a b]);
%
    % solution curve in phase space
    subplot(1,2,1);
    plot(ysol(:,1),ysol(:,2),'r-')
    hold on
%
    subplot(1,2,2);
    plot(t,ysol(:,1),'b-',t,ysol(:,2),'r-')
    legend('x','y')
    pause;
    end
end
%
disp('To be shown: solution curves with initial condition lies in [1,2]')
clf
subplot(1,2,1);
   quiver(x,y,u./leng,v./leng,0.4)
hold on
plot(double(x0),double(y0),'go')
for a=1.1:0.3:2
    for b = 1:0.4:2
    [t,ysol] = ode45(fode,[t0 tmax], [a b]);
%
    % solution curve in phase space
    subplot(1,2,1);
    plot(ysol(:,1),ysol(:,2),'r-')
%
    subplot(1,2,2);
    plot(t,ysol(:,1),'b-',t,ysol(:,2),'r-')
    legend('x','y')
    pause;
    end
end
hold off
