//*******************************************************
//
//  Resolucion de Ecuaciones Diferenciales Ordinarias
//
//*******************************************************
//  Problemas de condicion inicial
//*******************************************************
//
// Resolver la sgte. ecuacion diferencial
//
//  y'(t,y) = 30 - 5y    -- Ec. Diferencial
//  y = 1 para t = 0     -- Condicion Inicial
//  t = [0,10]           -- Tiempo de resolcuion
//*******************************************************
// En general y' = f(t,y)
//
// Met. de Euler: 
//      y(k+1) = y(k) + h * f(t,y)
//
// Met. de Euler Mejorado (Heun): 
//      y(k+1) = y(k) + h/2 *(k1 + k2)
//      k1 = f(t,y(k))
//      k2 = f(t+h,y(k)+h*f(t,y(k)))
//
//*******************************************************

clear

// Generando la solucion real
t = 0:0.001:10;
y_real = -5*exp(-5*t)+6;

// Condiciones iniciales
y_euler(1) = 1;
y_heun(1) = 1;

a= 0;    // Tiempo inicial
b= 10;   // Tiempo final

// Probar con diferentes valores de paso
h = 0.5;
//h = 0.5;
//h = 0.25;
//h = 0.1;

n = (b-a)/h;   // Numero de pasos

deff('yp=f(y)','yp=30-5*y');    // Definiendo la funcion 
t_sim(1) = a;

// Resolviendo usando el Met. de Euler
for i=1:n
    y_euler(i+1) = y_euler(i) + h*feval(y_euler(i),f);
    t_sim(i+1) = t_sim(i) + h;
end


// Resolviendo usando el Met. de Euler Mejorado
for i=1:n
    k1 = feval(y_heun(i),f);
    y_aux = y_heun(i) + h*k1;
    k2 = feval(y_aux,f);
    y_heun(i+1) = y_heun(i) + h/2*(k1+k2);
end


// Utilizando ahora la funcion de Scilab

deff('yp=f(t,y)','yp=30-5*y');    // Redefiniendo la funcion : convension Scilab f(t,y) (f debe ser funcion te t e y)
y_sci=ode("rk",1,a,t_sim,f)  

// Graficando los resultados
xbasc()   // Limpiando la ventana de graficos
plot2d(t,y_real)
plot2d(t_sim,[y_euler y_heun],[2,3],leg="Euler@Heun")
plot2d(t_sim,[y_euler y_heun y_sci' ],[2,3,4],leg="Euler@Heun@Scilab")