% This is code to draw the direction field and phase portrait
% for the system x'=-(x-y)(1-x-y); y'=x(2+y)

% first set the range by setting a,b,c,d
% to draw region [a,b]*[c,d] in the phase plane
a=-4;b=4;c=-4;d=4;

[x, y] = meshgrid(a:0.05:b, c:.05:d);

%% Notice that this part needs to be changed for different systems
xdot = -(x-y).*(1-x-y); %Note the use of .* and .^
ydot = x.*(2+y);
% Define the RHS function for numerical solutions later
RHS=@(t,z)[-(z(1)-z(2)).*(1-z(1)-z(2)); z(1).*(2+z(2))]; % z(1)=x;, z(2)=y


close all
plotdirectionfield=1; % if 0, turn off the direciton field
plottrajectory=1;

r=sqrt(xdot.^2+ydot.^2);
FigHandle = figure;
set(FigHandle, 'Position', [700, 500, 1249, 895]);
set(0,'DefaultTextFontSize',30,'DefaultAxesFontSize', 30)
axis([a b c d])
hold all;
% draw x-y axis
line([a b], [0 0],'color','k','linewidth',2) %x-axis
line( [0 0], [c d],'color','k','linewidth',2) %y-axis 
% draw critical points
plot(0,0,'*r', 'MarkerSize',12);
plot(0,1,'*r', 'MarkerSize',12);
plot(-2,-2,'*r', 'MarkerSize',12);
plot(3,-2,'*r', 'MarkerSize',12);

if (plotdirectionfield)
quiver(x,y,xdot./r, ydot./r); hold all
title('Normalized Direction Field')
end

%% to plot a few trajectories, we need to solve the system numerically
if (plottrajectory)
% draw a trajectory starting from (0.2,0.2) in the phase plane
[t z]=ode45(RHS,[0 20],[0.2; 0.2]); % compute a numerical solution from t=0 to t=20
plot(z(:,1),z(:,2),'-b'); hold all
% draw a trajectory starting from (1,1.2) in the phase plane
[t z]=ode45(RHS,[0 20],[1; 1.2]);
plot(z(:,1),z(:,2),'-b'); hold all
% draw a trajectory starting from (4,4) in the phase plane
[t z]=ode45(RHS,[0 20],[4; 4]);
plot(z(:,1),z(:,2),'-b'); hold all

% draw trajectories starting from N points on the x-axis on [-2,2] in the phase plane
N=10;
for i=1:N
    s(i)=-2+4/N*i; 
    [t z]=ode45(RHS,[0 20],[ s(i);0]);
    plot(z(:,1),z(:,2),'-b'); hold all
    
    [t z]=ode45(RHS,[0 -20],[ s(i);0]);
    plot(z(:,1),z(:,2),'-b'); hold all
end

% draw trajectories starting from N points on the y-axis on [-3,3] in the phase plane
N=10;
for i=1:N
    s(i)=-3+6/N*i; 
    [t z]=ode45(RHS,[0 20],[0;s(i)]);
    plot(z(:,1),z(:,2),'-b'); hold all
    
    [t z]=ode45(RHS,[0 -20],[0;s(i)]);
    plot(z(:,1),z(:,2),'-b'); hold all
end
end