%Competition Example from Class
%McDermott
%February 05, 2013
% Find equilibrium values, x=N1 and y=N2
syms x y
sys1=x*(4-x-y); %declaring 1st equation
sys2= y*(6-y-3*x); %declaring 2nd equation
[xc,yc]= solve(sys1, sys2, x,y) %solving the system of equations
disp('Critical Points:'); disp([xc,yc]); %displays the equilibrium points
%________________________________________________________________
%Phase diagram analysis
%x(1)=N1 x(2)=N2
warning off all %hides warning
f1=@(t,x)[x(1)*(4-x(1)-x(2));x(2)*(6-x(2)-3*x(1))]; %declaring the 2 differential equations
figure; hold on %keeps the graph to place future graphs on top of it (e.g. separatices, arrows)
%begin simulation
for a = -1.5:0.5:5 %a and b provide a starting place for the separatices- these are the red lines on the phase diagram
for b = -1.5:0.5:5 %these for loops help us draw the separatrices quicker
[t,xa] = ode45(f1, [0 20], [a b]); %solves the path over time 0-20 and starting points [a b]
plot(xa(:,1), xa(:,2),'r')
[t,xa] = ode45(f1, [0 -20], [a b]);
plot(xa(:,1), xa(:,2),'r')
end %ends loop for b
end %ends loop for a
hold on %retain previous graph while adding new graph
[X,Y] = meshgrid(0:0.2:4.5, 0:0.2:6.5); %creates a grid on which to place the arrows of motion or the meshgrid produces the x and y coordinates within the phase plane
%(0:0.2:4.5, 0:0.2:6.5) represent dividing the x grid from 0 to 4.5 every 0.2 spaces. Also, we divide the y axis by 0 to 6.5 every 0.2 spaces.
U = X.*(4-X-Y);
V = Y.*(6-Y-3*X);
L = sqrt((U/2).^2 + (V/4).^2); %helps draw the arrows
quiver(X,Y, U./L, V./L, 0.4); %draws the arrows of motion on the phase diagram
%quiver produces the direction/arrows for each pair of coordinates
%.4 give lengths of arrows
axis ([0 4.5 0 6.5]) %gives the axis length, x between 0 and 4.5 and y between 0 and 6.5
xlabel('N1') %label x axis
ylabel('N2') %label y axis
title 'Vector Field and Trajectories of N1 and N2' %Label the title
%___________________________________________________________________
%Solving ODE
hold off %provides a clean slate for new graphs
f1=@(t,x)[x(1)*(4-x(1)-x(2));x(2)*(6-x(2)-3*x(1))]; %declaring the differential equations
%let N1=x(1) and N2=x(2), Notice that the result is what we would expect
%from the phase diagram
ode45(f1, [0 20], [1 2], '-o') %solving and graphing the different equations over time 0-20 (i.e. [0 20]) using starting values of 1 and 2 (i.e. [1 2]).
%Also '-o' creates a thicker line when the ODE is solved.