//Euler modificado com regra do ponto médio 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 = -12*y
endfunction

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

//Definindo os parâmetros numéricos
dt=0.1; //Passo de tempo
n=20;     //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;

//Determinando o primeiro passo com Euler aprimorado
   yp(2) = y(1) + dt*(f(t(1),y(1)));
   t(2) = t(1) + dt;
   y(2) = y(1) + (dt/2)*((f(t(1),yp(2))) + f(t(1),y(1)));

//Resolvendo para os próximos passos de tempo
for i = 2:n
   t(i+1) = t(i) + dt;
   y(i+1) = y(i) + 2*dt*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]);