clear
clc

f = @(t,y) y^2 - y^3;

delta = 1e-5;
y0 = delta;
intervalo = [0 2/delta];
opts = odeset('AbsTol',1e-4, 'RelTol',1e-4);

[t,y] = ode23s(f, intervalo, y0,opts);

% comparar la cantidad de pasos que necesita cada solver!
size(y)
plot(t,y,'-o')


%% Usando Heun, necesitamos un monton de pasos!

cant_pasos = round(0.5/delta)
h = intervalo(2)/cant_pasos;

y = zeros(cant_pasos,1);
y(1) = y0;

k1 = 0;
k2 = 0;

for n = 1:cant_pasos
    tn = h*(n-1);
    k1 = h*f(tn, y(n));
    k2 = h*f(tn+h, y(n) + k1);
    y(n+1) = y(n) + (k1 + k2) /2;
end

t = h*(0:cant_pasos);
plot(t,y,'-o')