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

%Critical points:

global a b c d V k m n A;
b = 0.73;
k = 0.9;
c = 0.015;
d = 0.01;
timespan = 1000; %simulation time
initial_x_min = 0.4;
initial_x_max = 0.8;
numberOfPoint = 6; % number of starting points
yOverx = 1; %yOverx = intial value of x / intial value of y

syms x y;
vars = [x, y];
eqs = [x/(1+x)-b*x*y/(k+x), c*x*y/(x*y+1) - d*y];
%eqs = [a*x-b*x*y, c*x*y - d*y];
options = odeset('RelTol', 1e-6);
colors = 'rgbycmk';
figure, hold on
[xs,ys] = solve(eqs(1), eqs(2));
xc=double(xs)
yc=double(ys)
xn=numel(xs);
for point=1:1:xn
xc=double(xs(point));
yc=double(ys(point));
if ((xc>=0) && (yc>=0)) 
plot(xc,yc,'r.', 'MarkerSize',25);
text(xc+xc/100,yc,['[',num2str(xc),', ',num2str(yc),']']);
end;
end;

% xc=double(xs(2));
% yc=double(ys(2));
% Dlim=b*c/d*(1-k)*xc/(1+xc)/(k+xc)/(k+xc)
% if (xc>0) 
% plot(xc,yc,'r.', 'MarkerSize',25);
% text(xc+0.5,yc,['[',num2str(xc),', ',num2str(yc),']']);
% end;
% xc=double(xs(3));
% yc=double(ys(3));
% if (xc>0)
% plot(xc,yc,'r.', 'MarkerSize',25);
% text(xc+1,yc,['[',num2str(xc),', ',num2str(yc),']']);
% end;
counter=0;
for x0 = initial_x_min:(initial_x_max-initial_x_min)/(numberOfPoint-1):initial_x_max
    [t,X] = ode45(@lotkavolterra, [0, timespan], [x0;x0/yOverx],options);
    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.


figure, hold on
counter=0;
for x0 = initial_x_min:(initial_x_max-initial_x_min)/(numberOfPoint-1):initial_x_max
    [t, X] = ode45(@lotkavolterra, [0, timespan], [x0; x0/yOverx], options);
    subplot(2, 1, 1), hold on
    plot(t, X(:,1), colors(mod(counter,7)+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
