clear
clf

global alpha beta delta gamma h 
alpha = 1;
beta = 200;
delta = 300;
gamma = 1;


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

x0 = 400;
y0 = 100;

Y = [x0; y0]; 


% This part of the code is for your information only, regarding stable
% fixed points of the system
syms x y
vars = [x, y];
eqs = [x*(alpha-y/beta), -y*(gamma - x/delta) - h];
[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)
disp('Matrix of linearized system:')
subs(A, vars, [xc(2), yc(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_ode(t     , Y                );
    f2 = lotka_ode(t+dt/2, Y+f1 .*  (dt/2)  );
    f3 = lotka_ode(t+dt/2, Y+f2 .*  (dt/2)  );
    f4 = lotka_ode(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);

    k = animatedline('Color', 'b', 'linewidth', 2);
    j = animatedline('Color', 'm', 'linewidth', 2);
    axis([0, 30, 0, 500])
    title('Population of lemon sharks and mullet with infinite food')
    legend('Mullet','Lemon sharks')
    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));
drawnow
end
