// Ejemplito de metodos numericos para ecuaciones diferenciales ordinarias

function yu1=f1(t,x,y)
yu1 = y
endfunction

function yu2=f2(t,x,y)
k = 0.5;
l = 1.0;
yu2 = -k*y - l*sin(x);
endfunction


N  = 200; // Particion 
t0 = 0.0; // Tiempo inicial
t1 = 30.0; // Tiempo final
cix = 0.0; // Condicion inicial al tiempo t0
ciy = 3.0;

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

// Asigna la condicion inicial 
e1 = cix;
e2 = ciy;

// Solo para graficar
xgrid
X = zeros(N+1,1);
Y = zeros(N+1,1);
X(1) = cix;
Y(1) = ciy;

// Calculo de las soluciones
for i=1:N 
   tiempo = t+h // Solo para la visualizacion
 
   X(i+1)  = e1 + h * f1(t,e1,e2)  // Euler
   Y(i+1)  = e2 + h * f2(t,e1,e2)  

   e1 = X(i+1); // evolucion para graficar
   e2 = Y(i+1); 

   t = t0 + i * h; // Calcula en paso en el tiempo
end

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