%% Resolver el PVI y'(t) = a*y(t)*(1-y(t)), y(0) = 0.1. 
% La solucion exacta es y(t) = exp(t*a)./(exp(t*a)+1/y0-1)

clear 
clc
close

a = 1;
y0 = 0.1;
T = 8;

cant_pasos = 80;
h = T/cant_pasos;
y = zeros(1,cant_pasos+1);
y(1) = y0;
z  = y; % Euler h/atras

% Euler hacia adelante
for n = 1:cant_pasos
    y(n+1) = y(n) * (1 + a*(1-y(n))*h);
end


%% Euler hacia atras
for n = 1:cant_pasos
    disc = sqrt((1 - a*h)^2 + 4 * a *h *z(n));
    z(n+1) = (-(1 - a*h) + disc) / (2 * a * h);
end

%% Graficar
t = linspace(0,T);
plot(linspace(0,T,cant_pasos+1),y, 'ob', 'LineWidth', 2)
hold on
plot(linspace(0,T,cant_pasos+1),z, 'or', 'LineWidth', 2)
plot(t,exp(t*a)./(exp(t*a)+1/y0-1), '-k', 'LineWidth', 2)
grid on
xlabel('t', 'FontSize',14)
ylabel('y(t)', 'FontSize',14)

% y(end)- exp(t(end)*a)./(exp(t(end)*a)+1/y0-1)