%----------------------------%
% Lotka-Volterra equations   %
%----------------------------%

clear all
close all
clc

% dx/dt=a*x-b*x*y;
% dy/dt=-c*y+d*x*y;

a=1.5; b=1; c=3; d=1;

t0=0; tfin=20; N=2000;
t=linspace(t0,tfin,N);
dt=t(2)-t(1);

x(1)=10;
y(1)=5;

for i=1:N-1
    xm(i)=x(i) + (dt/2)*( a*x(i)-b*x(i)*y(i) );
    ym(i)=y(i) + (dt/2)*( -c*y(i)+d*x(i)*y(i) );
    
    x(i+1)=x(i) + dt*( a*xm(i)-b*xm(i)*ym(i) );
    y(i+1)=y(i) + dt*( -c*ym(i)+d*xm(i)*ym(i) );
end


figure(1)
hold on
plot(t,x,'r-','LineWidth',2)
plot(t,y,'b-','LineWidth',2)
legend('x(t)', 'y(t)',...
    'Location','NorthWest')
hold off

figure(2)
hold on
plot(x,y,'g-','LineWidth',2)
legend('phase space',...
    'Location','NorthEast')
hold off

%%
clear all
close all
clc

% dx/dt=a*x-b*x*y;
% dy/dt=-c*y+d*x*y;

a=1; b=1; c=1; d=1;

t0=0; tfin=40; N=5000;
t=linspace(t0,tfin,N);
dt=t(2)-t(1);

x(1)=10;
y(1)=5;

for i=1:N-1
    xm(i)=x(i) + (dt/2)*( a*x(i)-b*x(i)*y(i) );
    ym(i)=y(i) + (dt/2)*( -c*y(i)+d*x(i)*y(i) );
    
    x(i+1)=x(i) + dt*( a*xm(i)-b*xm(i)*ym(i) );
    y(i+1)=y(i) + dt*( -c*ym(i)+d*xm(i)*ym(i) );
end


figure(1)
hold on
plot(t,x,'r-','LineWidth',2)
plot(t,y,'b-','LineWidth',2)
legend('x(t)', 'y(t)',...
    'Location','NorthWest')
hold off

figure(2)
hold on
plot(x,y,'g-','LineWidth',2)
legend('phase space',...
    'Location','NorthEast')
hold off

%%
% using ode45

clear all
close all
clc

options = odeset('RelTol',1e-15,'AbsTol',[1e-11]);
[T,Y] = ode45(@equations, [0 40], [10 5]);


figure(1)
hold on
plot(T,Y(:,1),'r-','LineWidth',1)
plot(T,Y(:,2),'b-','LineWidth',1)
legend('x(t)', 'y(t)',...
    'Location','NorthWest')
hold off

figure(2)
hold on
plot(Y(:,1),Y(:,2),'g-','LineWidth',1)
legend('phase space',...
    'Location','NorthEast')
hold off
%%
% with an external function
clear all
close all
clc

% dx/dt=a*x-b*x*y;
% dy/dt=-c*y+d*x*y;

a=4; b=2; c=3; d=3;
% with these values, the stable center is at
xst=c/d 
yst=a/b

t0=0; tfin=10; N=2000;

figure(1)
hold on

x0=0.8;
y0=2.2;
[ x, y, t ] = RK_pray_pred( t0, tfin, N, x0, y0, a, b, c, d );

figure(1)
hold on
plot(t,x,'r-','LineWidth',2)
plot(t,y,'b-','LineWidth',2)
legend('x(t)', 'y(t)',...
    'Location','NorthWest')
hold off

figure(2)
hold on
plot(x,y,'g-','LineWidth',2)
legend('x0=0.8, y0=2.2',...
    'Location','NorthEast')
hold off

%%
% different initial conditions: look for the stable center
clear all
close all
clc

% dx/dt=a*x-b*x*y;
% dy/dt=-c*y+d*x*y;

a=4; b=2; c=3; d=3;
% with these values, the stable center is at
xst=c/d 
yst=a/b

t0=0; tfin=10; N=2000;

figure(1)
hold on

Nin=4
x0=linspace(0.1,3,Nin);
y0=linspace(0.1,3,Nin);

for i=1:Nin
    for j=1:Nin
        [ x, y, t ] = RK_pray_pred( t0, tfin, N, x0(i), y0(j), a, b, c, d );
        plot(x,y,'r-','LineWidth',1)
    end
end

f=@(x,y)(a*x-b*x.*y);
g=@(x,y)(-c*y+d*x.*y);

Nm=7;
xmin=0; xmax=3;
ymin=0; ymax=6;

xmesh=linspace(xmin,xmax,Nm);
ymesh=linspace(ymin,ymax,Nm);
[xp,yp]=meshgrid(xmesh,ymesh);

quiver(xp,yp,f(xp,yp),g(xp,yp),2.5)
%axis image

hold off

% x0=10;
% y0=10;
% 
% plot(x,y,'g-','LineWidth',2)
% 
% x0=0.1;
% y0=0.1;
% [ x, y, t ] = RK_pray_pred( t0, tfin, N, x0, y0, a, b, c, d );
% plot(x,y,'r-','LineWidth',2)
% 
% x0=1;
% y0=5;
% [ x, y, t ] = RK_pray_pred( t0, tfin, N, x0, y0, a, b, c, d );
% plot(x,y,'b-','LineWidth',2)
% 
% x0=0.8;
% y0=2.2;
% [ x, y, t ] = RK_pray_pred( t0, tfin, N, x0, y0, a, b, c, d );
% plot(x,y,'c-','LineWidth',2)
% 
% legend('x0=10, y0=10', 'x0=0.1, y0=0.1', 'x0=1, y0=5', 'x0=0.8, y0=2.2',...
%     'Location','NorthEast')

hold off
%%

%%
% different initial conditions: look for the stable center
clear all
close all
clc

% dx/dt=a*x-b*x*y;
% dy/dt=-c*y+d*x*y;

a=1; b=1; c=1; d=1;
% with these values, the stable center is at
xst=c/d 
yst=a/b

t0=0; tfin=10; N=2000;

figure(1)
hold on

Nin=4
x0=linspace(0.1,3,Nin);
y0=linspace(0.1,3,Nin);

plot(1,1,'sk','MarkerSize',4)
plot(0,0,'xk','MarkerSize',4)

f=@(x,y)(a*x-b*x.*y);
g=@(x,y)(-c*y+d*x.*y);

for i=1:Nin
    for j=1:Nin
        [ x, y, t ] = RK_pray_pred( t0, tfin, N, x0(i), y0(j), a, b, c, d );
        plot(x,y,'r-','LineWidth',1)
    end
end
legend('stable center','unstable critical point','orbits')

Nm=7;
xmin=0; xmax=3;
ymin=0; ymax=3;

xmesh=linspace(xmin,xmax,Nm);
ymesh=linspace(ymin,ymax,Nm);
[xp,yp]=meshgrid(xmesh,ymesh);

quiver(xp,yp,f(xp,yp),g(xp,yp),2)
%axis image

xlabel('x(prays)','FontSize',16)
ylabel('y(predators)','FontSize',16)
title('Prays and predators: phase space')
hold off

print -f -depsc 'Prays_Predators.eps'

%%
% Jaime problem
clear all
close all
clc

% dx/dt=f(x,y);   f(x,y)=60*x-4*x^2-3*x*y;
% dy/dt=g(x,y);   g(x,y)=42*y-2*y^2-3*x*y;


f=@(x,y) ( 60*x-4*x^2-3*x*y );
g=@(x,y) ( 42*y-2*y^2-3*x*y );

% critical points
% (0,0) unstable
% (6, 12) unstable
% (0,21) stable
% (15,0) stable

t0=0; tfin=5; N=1000;


Nin=5
x0=30;
y0=2;
[ x, y, t ] = RK_2coupledODE( t0, tfin, N, x0, y0, f, g );
figure(1)
hold on
plot(t,x,'r-','LineWidth',1)
plot(t,y,'b-','LineWidth',1)
hold off

%%
% Jaime problem: phase space
clear all
close all
clc
% dx/dt=f(x,y);   f(x,y)=60*x-4*x^2-3*x*y;
% dy/dt=g(x,y);   g(x,y)=42*y-2*y^2-3*x*y;


f=@(x,y) ( 60*x-4*x.^2-3*x.*y );
g=@(x,y) ( 42*y-2*y.^2-3*x.*y );

% critical points
% (0,0) unstable
% (6, 12) unstable
% (0,21) stable
% (15,0) stable

t0=0; tfin=10; N=2000;
figure(2)
hold on

Nin=5
x0=linspace(0,16,Nin);
y0=linspace(0,21,Nin);

for i=1:Nin
    for j=1:Nin
        [ x, y, t ] = RK_2coupledODE( t0, tfin, N, x0(i), y0(j), f, g );
        plot(x,y,'r-','LineWidth',1)
    end
end

% Add plot of f(x,y) and g(x,y) as vector fields

Nm=15;
xmin=0; xmax=14;
ymin=0; ymax=19;

xmesh=linspace(xmin,xmax,Nm);
ymesh=linspace(ymin,ymax,Nm);
[xp,yp]=meshgrid(xmesh,ymesh);

quiver(xp,yp,f(xp,yp),g(xp,yp),1.5)
%axis image

hold off

