
%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;

a = 0.4;
b = 0.01;
c = 0.01;
d = 0.4;



%m = 10;
%n = 30;
%V = 300;
% syms x y a b c d
% vars = [x, y]
% eqs = [x*(a-b*y), y*(-d+c*x)];
% [xc, yc] = solve(eqs(1), eqs(2));
% [xc, yc]
% A = jacobian(eqs, vars);
% disp('Matrix of linearized system:')
% subs(A, vars, [xc(1), yc(1)])
% disp('eigenvalues:')
% eig(ans)
% %eig(double(ans))
% disp('Matrix of linearized system:')
% subs(A, vars, [xc(2), yc(2)])
% disp('eigenvalues:')
% %eig(double(ans))
% eig(ans)

%phase plot

% eqs = @(t, x) ...
%       [x(1)*(-b*x(2)+V/(k+x(1))), x(2)*(-d+c*x(1))];
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-8);
colors = 'rgmcykb';
% figure, hold on
[xs,ys] = solve(eqs(1), eqs(2))
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',15);
% text(xc+0.5,yc,['[',num2str(xc),', ',num2str(yc),']']);
% text(xc+xc/10,yc,['[ d/c, a/b ]']);
end;
xc=double(xs(1));
yc=double(ys(1));
if (xc>0)
plot(xc,yc,'r.', 'MarkerSize',25);
text(xc+1,yc,['[',num2str(xc),', ',num2str(yc),']']);
end;
for x0 = 10:1:10
    [t,X] = ode45(@lotkavolterra, [0, 200], [x0;x0], options);
    plot(X(:,1), X(:,2), colors(mod(counter,7)+1),'LineWidth',2)
    
    counter=counter+1;
end, hold on

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 = 10:1:10
    [t, X] = ode45(@lotkavolterra, [0, 200], [x0; x0], options);
    subplot(2, 1, 1), hold on
    plot(t, X(:,1), colors(mod(counter,7)+1),'LineWidth',1.4)
    subplot(2, 1, 2), hold on
    plot(t, X(:, 2), colors(mod(counter,7)+1),'LineWidth',2)
    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
