%Modelo logístico ejercicio 8
% Euler
clear all
clf
%DATOS DEL PROBLEMA
t0=0;
tN=100;
y0=60;
h=input('introduce numero de paso: ');
M1=30;
M2=100;
r=0.04;
%Calculamos número de subintervalos
N=round((tN-t0)/h);
% Definimos la variable independiente
t=t0:h:tN;
%Ahora vamos a guardar los valores de la solución aproximada en el
%vector y
y=zeros(1,N+1); %euler
w=zeros(1,N+1); %Heun
z=zeros(1,N+1); %Runge-kutta
y(1)=y0;
w(1)=y0;
z(1)=y0;
for i=1:N
 y(i+1)=y(i)+h*(-r*y(i)*(1-(y(i)/M1))*(1-(y(i)/M2))); %euler
 k1=(-r*w(i)*(1-(w(i)/M1))*(1-(w(i)/M2))); %Heun
 k2=(-r*(w(i)+h*k1)*(1-((w(i)+h*k1)/M1))*(1-((w(i)+h*k1)/M2))); %Heun t(i)+h+w(i)+k1*h
 w(i+1)=w(i)+(h/2)*(k1+k2); %Heun
 % Runge-kutta
 A1=(-r*z(i)*(1-(z(i)/M1))*(1-(z(i)/M2)));
 A2=(-r*(z(i)+1/2*A1*h)*(1-((z(i)+1/2*A1*h)/M1))*(1-(z(i)+1/2*A1*h)/M2)); %z(i)+1/2*A1*h)*(1-(z(i)+1/2*A1*h))
 A3=(-r*(z(i)+1/2*A2*h)*(1-((z(i)+1/2*A2*h)/M1))*(1-(z(i)+1/2*A2*h)/M2));
 A4=(-r*(z(i)+A3*h)*(1-((z(i)+A3*h)/M1))*(1-(z(i)+A3*h)/M2));
 z(i+1)= z(i)+h/6*(A1+2*A2+2*A3+A4);
end
%sacamos tabla de resultados
[t',y',w',z']
hold on
%gráfico
plot(t,y,'r')
plot(t,w,'g')
plot(t,z)
legend('Euler','Heun','Runge-Kutta','Location','best'); % lo último es para que la leyenda
hold off