//****************************************************************
//
//  Resolucion de Sistema de Ecuaciones Diferenciales Ordinarias
//
//****************************************************************
//  Problemas de condicion inicial
//****************************************************************
//
// Resolver el sgte. sistema de ecuaciones diferenciales
//
//  x' = y               -- Ec. Diferencial
//  y' = -x + 6*cos(t)
//  x(0) = 2             -- Condicion Inicial
//  y(0) = 3
//  t = [0,10]           -- Tiempo de resolucion
//*******************************************************
// En general z' = f(t,z)   
// Como estamos resolviendo un sistema de ecuaciones
// z se trata de un vector z =[x y]
// en la primera componente se guarda el valor de x
// en la segunda componente se guarda el valor de y
//
// Met. de Euler: 
//      z(k+1) = z(k) + h * f(t,z)
//

//*******************************************************

clear


// Condiciones iniciales
z_inicial = [2 3];
z_euler = z_inicial;

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

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

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

getf('derivada_sistema.sci');    // Cargando la funcion donde se calcula la derivada

t_sim(1) = a;

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


// Utilizando ahora la funcion de Scilab

y_sci=ode("rk",z_inicial,a,t_sim,derivada_sistema)  


// Generando la solucion real (generalmente no conocida)(usada en este caso solo con fines de comparacion)
t = t_sim;
x_real = 2*cos(t)+3*sin(t)+3*t.*sin(t);
y_real = -2*sin(t)+3*cos(t)+3*sin(t)+3*t.*cos(t);


// Graficando los resultados
xset('window',1)
plot2d(t_sim,[x_real z_euler(:,1)],[2,3],leg="Real@Euler")
xtitle('Grafico x')

xset('window',2)
plot2d(t_sim,[y_real z_euler(:,2)],[2,3],leg="Real@Euler")
xtitle('Grafico y')