%% ------------------------Ej 9 --------------------------------------

% Parte a - sin evento
alpha = 0.01; % 0.1
y0 = [300 150];
F = @(t,y) [2*y(1) - alpha*y(1)*y(2); -y(2)+alpha*y(1)*y(2)];

[tiempos,poblaciones] = ode23(F,[0 1000],y0);

plot(poblaciones(:,1),poblaciones(:,2))
figure
plot(tiempos,poblaciones(:,1),'linewidth', 1), hold on
plot(tiempos,poblaciones(:,2),'linewidth', 1),
grid
legend('Conejos', 'Zorros')
hold off
%% 


%Parte a - con evento (mal)
alpha = 0.01;
y0 = [300 150];
F = @(t,y) [2*y(1) - alpha*y(1)*y(2); -y(2)+alpha*y(1)*y(2)];
opts = odeset('events',@parar9a);
[t,y,tfinal] = ode23(F,[0 10],y0,opts);

tfinal %Esto se calcula solo al usar ode23+evento.


plot(y(:,1),y(:,2))
figure
plot(t,y(:,1),'linewidth', 1), hold on
plot(t,y(:,2),'linewidth', 1),
grid
legend('Conejos', 'Zorros')
hold off






%%

%Parte a - con evento (bien)
alpha = 0.01;
y0 = [300 150];
F = @(t,y) [2*y(1) - alpha*y(1)*y(2); -y(2)+alpha*y(1)*y(2)];
opts = odeset('events',@parar_9a);
[t,y,tfinal] = ode23(F,[0 10],y0,opts);

tfinal


plot(y(:,1),y(:,2))
figure
plot(t,y(:,1),'linewidth', 1), hold on
plot(t,y(:,2),'linewidth', 1),
grid
legend('Conejos', 'Zorros')
hold off


%%

% Parte b - sin parar
alpha = 0.01;
y0 = [15 22];
F = @(t,y) [2*y(1) - alpha*y(1)*y(2); -y(2)+alpha*y(1)*y(2)];
[t,y] = ode23(F,[0 10],y0);

tfinal


plot(y(:,1),y(:,2))
figure
plot(t,y(:,1),'linewidth', 1), hold on
plot(t,y(:,2),'linewidth', 1),
grid
legend('Conejos', 'Zorros')
hold off

%% 

% Parte b - Parando

alpha = 0.01;
y0 = [15 22];
F = @(t,y) [2*y(1) - alpha*y(1)*y(2); -y(2)+alpha*y(1)*y(2)];
opts = odeset('events',@parar_9b);
%Cual va a ser la nueva regla de parada?
[t,y,tfinal] = ode23(F,[0 10],y0,opts);

tfinal


plot(y(:,1),y(:,2))
figure
plot(t,y(:,1),'linewidth', 1), hold on
plot(t,y(:,2),'linewidth', 1),
grid
legend('Conejos', 'Zorros')
hold off

%%
% Parte c

alpha = 0.01;
y0 = [102 198];
F = @(t,y) [2*y(1) - alpha*y(1)*y(2); -y(2)+alpha*y(1)*y(2)];
[t,y] = ode23(F,[0 10000],y0);

plot(y(:,1),y(:,2))
figure
plot(t,y(:,1),'linewidth', 1), hold on
plot(t,y(:,2),'linewidth', 1),
grid
legend('Conejos', 'Zorros')
hold off


%%

%parte c- parando?
alpha = 0.01;
y0 = [102 198];
F = @(t,y) [2*y(1) - alpha*y(1)*y(2); -y(2)+alpha*y(1)*y(2)];
opts = odeset('events',@parar_9c);
[t,y,tfinal] = ode23(F,[0 10],y0,opts);

tfinal


plot(y(:,1),y(:,2))
figure
plot(t,y(:,1),'linewidth', 1), hold on
plot(t,y(:,2),'linewidth', 1),
grid
legend('Conejos', 'Zorros')
hold off

%% -------------------------- Ej 3 ------------------------------------------

%Parte b)1)

% Definimos la ecuación dy/dt = f(t, y)
f = @(t, y) -y^2; % Funcion

% t0, tf, y0 y h
t0 = 0;
tf = 1;
y0 = 2;
h = 0.1;

% Llamada a euler
[t, y] = eulerE(f, t0, tf, y0, h);

% Plot 
plot(t, y, '-o');
xlabel('t');
ylabel('y');
title('Euler directo para f');
legend('Euler directo');
hold on;


%% Con Euler implicito

% Parte b)2)

%  dy/dx = f(t, y)
f = @(t, y) -y^2; % Funcion
df_dy = @(t, y) -2*y; % Derivada (para newton)

% Cond iniciales
t0 = 0;
tf = 1;
y0 = 2;
h = 0.1;

% Euler Implicito
[t, y] = EulerI(f, df_dy, t0, tf, y0, h)

% Plot 
plot(t, y, '-o');
xlabel('t');
ylabel('y');
title('Euler Implicito para f');
legend('Euler Implicito');
hold on;


%% Con metodo del trapecio

% Ecuacion dy/dx = f(t, y)
f = @(t, y) -y^2; % funcion f
df_dy = @(t, y) -2*y; % Derivada parcial para pasarsela a la funcion de trapecio
% Condiciones iniciales
t0 = 0;
tf = 1;
y0 = 2;
h = 0.1;

% Llamada a la función
[t, y] = trapecio(f, df_dy, t0, tf, y0, h);

% Plot
plot(t, y, '-o');
xlabel('t');
ylabel('y');
title('Método del trapecio');
legend('Trapecio');


%% -------------------------  Ejercicio 4 ---------------

f = @(t, y) -2*t*(1+t)*y.^2 - y./(1+t); % la ode dy/dt=f
df_dy = @(t, y) -4*t*(1+t)*y - 1./(1+t); % derivada parcial (con resp de y)

% Condiciones iniciales
t0 = 0;
tf = 1;
y0 = 1;
h = 0.1;

% Llamada a la función
[t, y] = trapecio(f, df_dy, t0, tf, y0, h);
error = y(11)-1/4 % y(1)=1/4.
% Plot
plot(t, y, '-o');
xlabel('t');
ylabel('y');
title('Método del trapecio');
legend('Trapecio');



%% -------------------------------Functions----------------------------


%% Parar los solvers 


function [val,isterm,dir] = parar9a(t,y)
    val = y(2) - 150; % y(1)-300;
    isterm = 1;
    dir = [];
end



function [val,isterm,dir] = parar_9a(t,y) %busca el primer momento donde y(1) = 300 y viene subiendo la funcion.
    val = y(1)-300;
    isterm = 1;
    dir = [1]; % si quiero que decrezca dir = [-1];
end

function [val,isterm,dir] = parar_9b(t,y)
    val = y(2)-22;
    isterm = 1;
    dir = [-1];
end

function [val,isterm,dir] = parar_9c(t,y)
    val = y(2)-22;
    isterm = 1;
    dir = [-1];
end

%------------------- Euler Directo ----------------------

function [t, y] = eulerE(f, t0, tf, y0, h)
    % Numero de divisiones
    n = floor((tf - t0) / h);
    
    % Creo arrays para guardar las sols
    t = zeros(1, n+1);
    y = zeros(1, n+1);
    
    % Cond iniciales
    t(1) = t0;
    y(1) = y0;
    
    % Implementacion de la iteracion
    for i = 1:n
        t(i+1) = t(i) + h;
        y(i+1) = y(i) + h * f(t(i), y(i));
    end
end

% ------------------ Euler Implicito -----------------

function [t, y] = EulerI(f, df_dy, t0, tf, y0, h)
    % numero de divisiones
    n = floor((tf - t0) / h);
    
    % Creo arrays para guardar las sols
    t = zeros(1, n+1);
    y = zeros(1, n+1);
    
    % Condiciones iniciales
    t(1) = t0;
    y(1) = y0;
    
    % Euler hacia atras
    for i = 1:n
        t(i+1) = t(i) + h;
        % adivino y(i+1)
        y_guess = y(i);
        
        % Newton para resolver la ecuación implícita
        tol = 1e-6; % Tolerancia
        max_iter = 10; % Numero max de iter
        for k = 1:max_iter
            F = y_guess - y(i) - h * f(t(i+1), y_guess);
            J = 1 - h * df_dy(t(i+1), y_guess); % Jacobiano (der resp y)
            
            if abs(J) < tol
                error('Jacobiano chico, puede haber problemas.');
            end
            
            y_new = y_guess - F / J; % Actualizar la variable
            if abs(y_new - y_guess) < tol
                y_guess = y_new;
                break;
            end
            y_guess = y_new;
        end
        
        y(i+1) = y_guess;
    end
end

% ------------ Trapecio  -----------------


function [t, y] = trapecio(f, df_dy, t0, tf, y0, h)
    % Numero de divisiones
    n = floor((tf - t0) / h);
    
    % Creo arrays para guardar las sols
    t = zeros(1, n+1);
    y = zeros(1, n+1);
    
    % Cond Iniciales
    t(1) = t0;
    y(1) = y0;
    
    % Implemento el trapecio
    for i = 1:n
        t(i+1) = t(i) + h;
        
        % adivino y(i+1)
        y_next = y(i);
        
        % Newton para resolver la ecuación implícita
        tol = 1e-6; % Tolerancia
        max_iter = 3; % Numero max de iter
        for k = 1:max_iter
            F = y_next - y(i) - (h/2) * (f(t(i), y(i)) + f(t(i+1), y_next));
            J = 1 - (h/2) * df_dy(t(i+1), y_next); % Jacobiano            
            if abs(J) < tol
                error('Jacobiano chico, puede haber problemas.');
            end
            
            y_new = y_next - F / J; % Paso de actualizacion
            if abs(y_new - y_next) < tol
                y_next = y_new;
                break;
            end
            y_next = y_new;
        end
        
        y(i+1) = y_next;
    end
end