clc; close; clear;

f=@(t,y) 3*y*(1-y)-y.^2/(1+y.^2);
ft=@(t,y) 0;
fy=@(t,y) 3-6*y-2*y./(1+y.^2).^2;

t0=0; tf=10;
haux=[0.25 0.5 1 2]; %tamanho da partição
y0=0.1;

for j=1:4
subplot(2,2,j)
    h=haux(j);
    n=(tf-t0)/h;
    t=t0:h:tf;
%Euler  
y=y0*ones(size(t));
    for i=1:n
       y(i+1)=y(i)+h*f(t(i),y(i));
    end
plot(t,y,'-*','LineWidth',2)
hold on
%Taylor2
y=y0*ones(size(t));
for i=1:n
    y(i+1)=y(i)+h*f(t(i),y(i))+((h^2)/2)*(ft(t(i),y(i))...
                       +fy(t(i),y(i))*f(t(i),y(i)));
end
plot(t,y,'-*','LineWidth',2)
%RK2-Heun
y=y0*ones(size(t));
for i=1:n
   K1=f(t(i),y(i));
   yp=y(i)+K1*h;
   K2=f(t(i)+h,yp);
   y(i+1)=y(i)+h*(K1+K2)/2;
end
plot(t,y,'-*','LineWidth',2)
%RK4
y=y0*ones(size(t));
for i=1:n
   K1=f(t(i),y(i));
   K2=f(t(i)+h/2,y(i)+K1*h/2);
   K3=f(t(i)+h/2,y(i)+K2*h/2);
   K4=f(t(i)+h,y(i)+K3*h);
   y(i+1)=y(i)+h*(K1+2*K2+2*K3+K4)/6;
end
plot(t,y,'-*','LineWidth',2)

legend('Euler','Taylor 2','Heun (RK2)','RK4')
ylim([0 1.2])
grid on
title(['h=',num2str(h)])
end

