Sistemas de Ecuaciones diferenciales con Octave

Este es el post que realmente quería escribir.  Pueden visitar el anterior post para ver como se resuelve una ecuación diferencial con Octave, ahora les explicare como hacerlo para resolver un sistema de ecuaciones diferenciales por método numérico, utilizando para ello un ejemplo con el modelo SIR (para epidemias) propuesto por W. O. Kermack y A. G. McKendrick, donde

• S(t) representa a los individuos susceptibles, es decir, aquellos que no han enfermado anteriormente y por lo tanto pueden resultar infectados al entrar en contacto con la enfermedad.
• I(t) representa a los individuos infectados y por lo tanto en condiciones de transmitir la enfermedad a los del grupo S.
• R(t) representa a los individuos recobrados de la enfermedad, y que ya no están en condiciones ni de enfermar nuevamente ni de transmitir la enfermedad a otros.

Las ecuaciones diferenciales son:

dS = -(Constante_beta)SI    dI = (Constante_beta)SI - (Constante_gama)I
dt                                         dt

          dR = (Constante_gama)I
          dt

donde:  (Constante_beta) = 0.2  Constante_gama = 0.7 y las condiciones iniciales son:  S(0) = 19    I(0) = 1 R(0) = 0.

Ahora para resolverlo con lsode se realiza lo siguiente:
voy a colocar el código y luego lo explicare : )

octave:1> function sis = f(x,t); b=0.2;g=0.7;
> sis(1) = -b*x(1)*x(2);
> sis(2) = b*x(1)*x(2)-g*x(2);
> sis(3) = g*x(2);
>endfunction
octave:2> t = (0:0.1:5)';
octave:3> x = lsode("f",[19,1,0],t);
octave:4> plot(t,x)

obtenemos esto

Ahora, explicaré como es q funciona.  Cada ecuación diferencial es una ecuación sis(n).  cada resolución de la ecuación, se representa como x(n), así se va colocando cada ecuación diferencial del sistema, como por ejemplo para los individuos susceptibles dS/dt se escribio la siguiente ecuación:
sis(1) = -b*x(1)*x(2); para escribir: dS/dt = b*S*I, así x(1) representa S y x(2) a I. El otro cambio es para colocar las condiciones iniciales en lsode, donde se coloca en una lista.    Notese que se grafican las tres ecuaciones.  Para colocarle nombre y color para cada gráfica ya son comandos que se conocen en GNUPlot.

Ojalá les sirva estimados lectores : ) porque a mi, si me servira xD : )