//Euler aprimorado para resolver o PVI y'(t) = f(t,y), y(t0)=y0)

clear //Limpar a memória

//Definindo f(t,y))
function a=f(t,y)
    a =(1+2*t)*sqrt(y)
endfunction

//Definindo as condições iniciais
t0=1;
y0=1;

//Definindo os parâmetros numéricos
dt=0.25; //Passo de tempo
n=4;     //número de passos avaliados


//Alocando as condições iniciais na posição inicial (1) dos vetores t(i) e y(i) 
t(1)=t0;
y(1)=y0;

//Resolvendo para os próximos passos de tempo
for i = 1:n
   yp(i+1) = y(i) + dt*(f(t(i),y(i)));
   t(i+1) = t(i) + dt;
   y(i+1) = y(i) + (dt/2)*((f(t(i+1),yp(i+1))) + f(t(i),y(i)));
end

c=t0:dt:t0+dt*n
//Solução exata para comparação, quando possível:
exat = ode(y0,t0,c,f)

//Plotando os valores obtidos e a solução exata para comparação, quando possível
plot(t,y,c,exat);
legend(['Numérico','Exata'],[,-1]);