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

%Critical points:
clear all;
global a b c d V k m n A e;

b = 0.45;
k = 0.35;
c = 1;
d = 0.5;
e = 0.5;


alpha=b*c/d/e+k+(e-1)/e
beta=b*(c/d-1)/e+k*(e-1)/e
delta=b/e

timespan = 500;    % simulation time
initial_x_min = 4;
initial_x_max = 16;
numberOfPoint = 6;  % number of starting points
yOverx = 50;         % yOverx = intial value of x / intial value of y

syms x y;
vars = [x, y];
eqs = [x/(1+x)-b*x*y/(k+x)-e*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)
% Dlim=c/d*xc/yc/(k+xc)*((1-k)/(1+xc)/(1+xc)-e)
xn=numel(xs);
for point=1:1:xn
xc=double(xs(point));
yc=double(ys(point));
% temp = (1-k)/(1+xc)/(1+xc)
% Dlim=c/d*xc/yc/(k+xc)*((1-k)/(1+xc)/(1+xc)-e)
if ((xc>=0) && (yc>=0)) 
plot(xc,yc,'r.', 'MarkerSize',25);
text(xc+xc/100,yc,['[',num2str(xc),', ',num2str(yc),']']);
end;
end;

counter=0;
for x0 = initial_x_min:(initial_x_max-initial_x_min)/(numberOfPoint-1):initial_x_max
    [t,X] = ode45(@lotkavolterra2, [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'

figure, hold on
counter=0;
for x0 = initial_x_min:(initial_x_max-initial_x_min)/(numberOfPoint-1):initial_x_max
    [t, X] = ode45(@lotkavolterra2, [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