// Ejemplito de metodos numericos para ecuaciones diferenciales ordinarias

function yy=f(t,y)
yy = y - t * t + 1;
endfunction


N  = 10; // Particion 
t0 = 0.0; // Tiempo inicial
t1 = 2.0; // Tiempo final
ci = 0.5; // Condicion inicial al tiempo t0


h = (t1 - t0) / N; // incremento en el tiempo
t = t0; // Asigna el tiempo inicial

// Asigna la condicion inicial 
e  = ci; // Para Euler
pm = ci; // Para Punto medio
em = ci; // Para Euler modificado
mh = ci; // Para Heun

// Solo para graficar
xgrid
X = zeros(11,1);
Y = zeros(11,1);
Ye = zeros(11,1);
X(1) = t;
Y(1) = ci;
Ye(1) = ci;

// Calculo de las soluciones
for i=1:N 
   tiempo = t+h // Solo para la visualizacion
   X(i+1) = tiempo; // tiempo para graficar

   e  = e + h * f(t,e)  // Euler
   pm = pm + h * f(t + (h/2.0), pm + (h/2.0)* f(t,pm)) // Punto medio 
   em = em + (h/2.0)*(f(t,em)+ f(t + h, em + h * f(t,em))) // Euler modificado
   mh = mh + (h/4.0)*(f(t,mh)+ 3.0*f(t+(2./3.0)*h,mh+(2.0/3.0)*h*f(t,mh))) // Heun
   
   Y(i+1) = e; // evolucion para graficar

   t = t0 + i * h; // Calcula en paso en el tiempo
   exacta = (t+1)^2 - 0.5*exp(t) // Calcula la solucion exacta al tiempo t
   Ye(i+1) = exacta; // evolucion exacta para graficar
end

// Grafica la solucion
plot2d(X,Ye,5) // Exacta
plot2d(X,Y,3)  // Calculada