clear, clc

r1 = 0.05;
r2 = 0.08;
c = 3000;
c2 = 15000;
K1 = 150000;
K2 = 400000;
a = 10^-8;

% Find equilibrium points %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
syms x y

blue = r1*x*((x - c)/(x + c))*(1 - x/K1) - a*x*y;
fin = r2*y*((y - c2)/(y + c2))*(1 - y/K2) - a*x*y;

[xstar, ystar] = vpasolve(blue,fin,x,y);

% Make vector field %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[X, Y] = meshgrid(0:20000:150000, 0:20000:400000); 
% Note: We don't want to make the step size too small, since it could make 
% these matrices prohibitively large for Matlab to do operations

f1 = r1.*X.*((X - c)./(X + c)).*(1 - X./K1) - a.*X.*Y; 
f2 = r2.*Y.*((Y - c2)./(Y + c2)).*(1 - Y./K2) - a.*X.*Y;
% .* and ./ indicate "element-wise" operations, instead of multiplying the 
% whole matrix we just want to multiply each element with each element

% Make nullclines %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[Xn, Yn] = meshgrid(0:1000:150000, 0:1000:400000); % Need a finer mesh for the nullclines 
f1n = r1.*Xn.*((Xn - c)./(Xn + c)).*(1 - Xn./K1) - a.*Xn.*Yn; 
f2n = r2.*Yn.*((Yn - c2)./(Yn + c2)).*(1 - Yn./K2) - a.*Xn.*Yn;

% Plot %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure
hold on
quiver(X,Y,f1,f2)
contour(Xn,Yn,f1n,[0,0],'k')
contour(Xn,Yn,f2n,[0,0],'k')