clear
clf

global k1 M1 k2 M2 c c2 gamma delta zeta
k1 = .01;
M1 = 100;
k2 = .04;
M2 = 80;
c = .01;
c2 = .03;
zeta = .02;
delta = 100;
gamma = 2;


t    =  0.0;
tmax = 15;
dt   = (1/9);
n    = tmax./dt; % number of steps

x0 = 90;
y0 = 100;
z0 = 80;

Y = [x0; y0; z0]; 


% This part of the code is for your information only, regarding stable
% fixed points of the system
syms x y z
vars = [x, y, z];
eqs = [x.*(k1.*(M1-x) - c.*y) , -y.*(gamma - (x./delta) - zeta.*z), z.*(k2.*(M2 - z) - c2.*y)];
[xc, yc, zc] = solve(eqs(1), eqs(2), eqs(3));
[xc, yc, zc]
A = jacobian(eqs, vars);
disp('Matrix of linearized system:')
subs(A, vars, [xc(1), yc(1), zc(1)])
disp('eigenvalues:')
eig(ans)
disp('Matrix of linearized system:')
subs(A, vars, [xc(2), yc(2), zc(2)])
disp('eigenvalues:')
eig(double(ans))

tt = zeros(n+1,1);
Yt = zeros(n+1,size(Y,1));
tt(1)   = t;
Yt(1,:) = Y';

for i=1:n

    f1 = lotka_ode2(t     , Y                );
    f2 = lotka_ode2(t+dt/2, Y+f1 .*  (dt/2)  );
    f3 = lotka_ode2(t+dt/2, Y+f2 .*  (dt/2)  );
    f4 = lotka_ode2(t+dt  , Y+f3 .*   dt     );
   
    Y = Y + (f1 + f2.*2 + f3.*2 + f4) .* (dt/6);
    

    t = t+dt;
    tt(i+1)   = t;  % this records the time at each step
    Yt(i+1,:) = Y'; % this records the vector Y at each step --> Yt(:,1) has all x coordinates
    
end

xx = Yt(:,1);
yy = Yt(:,2);
zz = Yt(:,3);

    k = animatedline('Color', 'b', 'linewidth', 2);
    j = animatedline('Color', 'm', 'linewidth', 2);
    l = animatedline('Color', 'g', 'linewidth', 2);
    axis([0, 15, 0, 120])
    title('Population of lemon sharks, mullet, and jacks')
    legend('Mullet','Lemon sharks', 'Jacks')
    xlabel('Time (years)') % x-axis label

for i = 1:length(xx)
hold on
    addpoints(k, tt(i),  xx(i));
    addpoints(j, tt(i),  yy(i));
    addpoints(l, tt(i),  zz(i));
drawnow
end
