Octave allows solving the equations of a kinetic system, in a very simple and efficient way. We will apply Octave to two complex kinetic systems: consecutive reactions and oscillating reactions (Lotka mechanism)

**a) Consecutive reactions:**

$A\stackrel{k_1}{\rightarrow} B \stackrel{k_2}{\rightarrow} C$, we solve the following differential equations with Octave: \begin{eqnarray} \frac{d[A]}{dt}& = & -k_1[A]\\ \frac{d[B]}{dt}& = & k_1[A]-k_2[B]\\ \frac{d[C]}{dt}& = & k_2[B ] \end{eqnarray}

We create a function that gives the velocity of the three components as a function of concentrations:

octave> function xdot = f(x,t) # x represents the concentrations of A, B and C, t is the time at a given instant and xdot are the derivatives of [A], [B] and [C] with respect to the time.

> global k # k is a vector representing the kinetic constants $k_1$ and $k_2$ is defined globally to be visible in the function above.

> xdot(1) = -k(1)*x(1) ; # define the first equation

> xdot(2) = k(1)*x(1) - k(2)*x(2); # define the second equation

> xdot(3) = k(2)*x(2); # define the third equation.

> end function

Now we go on to indicate the values of rate constants and initial concentrations:

octave> global k; k(1) = 1.0; k(2) = 1.0; # define k(1) and k(2)

octave> conc0(1) = 5.0; # define the initial concentration of A

octave> conc0(2) = 0.0; # define the initial concentration of B

octave> conc0(3) = 0.0; # defines the initial concentration of C

Next, we indicate the times between which we are going to carry out the integration and the number of stages that octave must carry out during the integration.

octave> t_initial = 0.0;

octave> t_final = 10.0;

octave> n_steps = 200;

octave> t = linspace(t_start, t_end, n_steps); **#** Integration of the differential equations is done by calling the lsode routine included in octave.

octave> conc = lsode("f", conc0, t); **#** The output is the matrix conc, whose columns are concentrations of the three components and time. However, a plot of concentrations versus times is more visual. We write the code that makes the graph:

octave> gset xlabel "time";

octave> gset ylabel "concentrations";

octave> plot(t, conc);

**Lotka mechanism**

$A+X\stackrel{k_1}{\rightarrow} 2X; \;\;\;\; X+Y\stackrel{k_2}{\rightarrow} 2Y; \;\;\;\; Y\stackrel{k_3}{\rightarrow}B$

The equations that give us the variations in time of the different reactants are: \begin{eqnarray} \frac{d[A]}{dt} & = & 0.0\\ \frac{d[X]}{dt} & = & k_1[A][X]-k_2[X][Y]\\ \frac{d[Y]}{dt} & = & k_2[X][Y]-k_3[Y]\\ \frac{ d[B]}{dt} & = & k_3[Y] \end{eqnarray}

We solve these equations numerically using octave.

octave> function xdot = f(x,t) # x represents the concentrations of A, X , Y and B, t is the time at a given instant and xdot are the derivatives of [A], [X], [Y] and [B] with respect to time.

> global k

> xdot(1) = 0.0; # define the first function

> xdot(2) = k(1)*x(1) *x(2) - k(2)*x(2)*x(3); # define the second function.

> xdot(3) = k(2)*x(2)*x(3) - k(3)*x(3); # define the third equation.

> xdot(4) = k(3)*x(3); # define the fourth equation.

> end function

octave> global k; k(1) = 0.05; k(2) = 0.1; k(3) =0.05; #$ define k(1) and k(2)

octave> conc0(1) = 5.0; # define the initial concentration of A

octave> conc0(2) = 5.0e-4; # define the initial concentration of X

octave> conc0(3) = 1.0e-5; #define the initial concentration of Y

octave> conc0(3) = 0.0; # define the initial concentration of B

octave> t_initial = 0.0;

octave> t_final = 650.0;

octave> n_steps = 10000;

octave> t = linspace(t_start, t_end, n_steps);

octave> conc = lsode("f", conc0, t);

octave> plot(t, conc);