
%d prey/dt = a *prey -  b * prey * predator
% d predator/dt = c * prey * predator - d * predator

%Critical points:
%b>c+cd+bd
global a b c d; 
a = 1;
b = 1;
c = 0.5;
d = 0.3;

syms x y;
vars = [x, y];
eqs = [x*(1-x)-a*x*y/(d+x), b*x*y/(x+d) - c*y];

[xs,ys] = solve(eqs(1), eqs(2))
xc=double(xs(2));
yc=double(ys(2));

options = odeset('RelTol', 1e-6);
colors = 'rgbycmk';
figure, hold on
[xs,ys] = solve(eqs(1), eqs(2))
xc=double(xs(2));
yc=double(ys(2));
plot(xc,yc,'r.', 'MarkerSize',25);
text(xc+1,yc,['[',num2str(xc),', ',num2str(yc),']']);
xc=double(xs(3));
yc=double(ys(3));
plot(xc,yc,'r.', 'MarkerSize',25);
text(xc+1,yc,['[',num2str(xc),', ',num2str(yc),']']);
counter=0;
for x0 = 0.2:0.4:1.8
    [t,X] = ode45(@ode, [0, 250], [x0;x0]);
    plot(X(:,1), X(:,2), colors(mod(counter,7)+1))
    counter=counter+1;
end, hold off

title 'population of predator against prey'
xlabel 'x = prey'
ylabel 'y = predator'

%Plot of Populations vs. Time

%We color-code the plots so you can see which ones go together.

counter=0;
figure, hold on
for x0 =0.2:0.4:1.8
    [t, X] = ode45(@ode, [0, 250], [x0; x0], options);
    subplot(2, 1, 1), hold on
    plot(t, X(:,1), colors(mod(counter,7)+1))
    counter=counter+1;
    subplot(2, 1, 2), hold on
    plot(t, X(:, 2), colors(mod(counter,7)+1))
    counter=counter+1;
    hold on
end
subplot(2, 1, 1)
title 'population of prey or predator against time'
xlabel t
ylabel 'x = prey'
subplot(2, 1, 2)
xlabel t
ylabel 'y = predators'
hold off
