//RK4 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
    k1 = f(t(i),y(i));
    t(i+1) = t(i) + dt;
    k2 = f((t(i)+dt/2),(y(i)+dt*k1/2));
    k3 = f((t(i)+dt/2),(y(i)+dt*k2/2));
    k4 = f(t(i+1),y(i)+dt*k3);
    y(i+1) = y(i) + (dt/6)*(k1 + 2*k2 + 2*k3 + k4)
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]);