% Sistema de 2 Ecuaciones Diferenciales
% Modelo "Lotka-Volterra" -> Crec. Logístico Presa
 
clc
close all
clear all

disp('Sistema de 2 Ecuaciones Diferenciales')
disp('Modelo "Lotka-Volterra" -> Crec. Logístico Presa')
 
disp('...................................................................')
 
disp('Ecuaciones: ')
 
disp('Presa -> dx/dt = a*x*(K-x)/K - b*x*y') 
disp('Depredador -> dy/dt = -c*y + d*x*y')

disp('...................................................................')
 
%disp('Valores por defecto:')
%disp('Valor inicial de Presas (x) x0: 3')
%disp('Valor inicial de Depredadores (y) y0: 1')
%disp('Valor de K: 20')
%disp('Tasa de nacimiento de la Presa, a: 0.4')  
%disp('Tasa de mortalidad del Depredador, c: 0.37') 
%disp('Susceptibilidad de la Presa a la Depredación, b: 0.3')  
%disp('Habilidad predatoria del Depredador, d: 0.05') 
%disp('tiempo final, tf: 100')
 
%disp('-------------------------------------------------------------------')

z=zeros(1,2);

z(1)=3; 
z(2)=1; 
K=20;
a=1;  
c=0.37; 
b=0.3;  
d=0.05; 
tf=100;

%z(1)=input('Valor inicial de x (x0): ');
%z(2)=input('valor inicial de y (y0): ');
%K=input('K: ');
%a=input('a: ');
%c=input('c: ');
%b=input('b: ');
%d=input('d: ');
%tf=input('tiempo final, tf: ');

%disp('...................................................................')

% z es una Matriz de 2 Columnas
% z(:,1) (todos los elementos de la 1era. Columna) representará los sucesivos valores de la función variable "x"
% z(:,2) (todos los elementos de la 2da.  Columna) representará los sucesivos valores de la función variable "y"

% Se definen las funciones "f(t,x,y)" y "g(t,x,y)" -> Ecuaciones Diferenciales "dx/dt" y "dy/dt" 
fg=@(t,z) [a*z(1)*((K-z(1))/K)-b*z(1)*z(2);-c*z(2)+d*z(1)*z(2)]; % z(1) es "x", z(2) es "y"

tspan=[0 tf]; 
[t,z]=ode23s(fg,tspan,z);

%**************************************************************************

% "x" vs. "t"
%figure
%plot(t,z(:,1),'r','LineWidth',2)
%grid
%xlabel('t (años)')
%ylabel('x (miles)');
%title('x vs. t')

%**************************************************************************

% "y" vs. "t"
%figure
%plot(t,z(:,2),'b','LineWidth',2)
%grid
%xlabel('t (años)')
%ylabel('y (miles)');
%title('y vs. t')

%**************************************************************************

% "x" e "y" vs. "t"
figure
plot(t,z(:,1),'r',t,z(:,2),'b','LineWidth',2)
grid
xlabel('t (años)')
ylabel('x,y (miles)');
legend('x','y');
title('x e y vs. t')

%**************************************************************************

% y vs. x  (Diagrama de Fases)
figure
plot(z(:,1),z(:,2),'Color',[0 0.5 0],'LineWidth',2)
grid
xlabel('x (miles)');
ylabel('y (miles)');
title('y vs. x (Diagrama de Fases)')

hold on
#line es una funcion util que esencialmente grafica la diagonal de un rectangulo, asi que se pone primero elrango en x y despues el rango en y
#en este caso directamente ponemos los valores que ya conocemos para nuestras nulclinas
#Observar que podrian ser mas largas, pero para visualizar asi esta bien
line([0 K],[a/b 0],'color','red'); 
line([c/d c/d],[0 a/b] ,'color','black');

line([0 0],[0 a/b],'color','red'); 
line([c/d c/d],[0 0] ,'color','black');

scatter(z(1,1),z(1,2),'filled') %marca el punto inicial, para que quede claro en qué sentido progresa el sistema
legend('(N,P)','dN/dt=0','dP/dt=0')



#Otras opciones para el diagrama de fase, descomentar para usar
#% y vs. x  (Diagrama de Fases 2)
#figure
#plot(z(:,1),z(:,2),'o','MarkerSize',2,'MarkerEdgeColor','k','MarkerFaceColor','r')
#%plot(z(:,1),z(:,2),'g','LineWidth',2)
#grid
#xlabel('x (miles)');
#ylabel('y (miles)');
#title('y vs. x (Diagrama de Fases)')