Revisión del 00:53 6 mar 2015

1 Introducción: objetivos y metodología

El objetivo de este trabajo es estudiar las concentraciones de los reactivos y productos de una reacción química a lo largo de la evolución en el tiempo de dicha reacción. Para ello modelizaremos el proceso empleando ecuaciones diferenciales y teniendo en cuenta la Ley de Acción de Masas y el Principio de Conservación de la materia.

Al ser las ecuaciones y sistemas que se nos presentan de difícil resolución analítica, utilizaremos el programa informático Matlab para resolverlos mediante métodos numéricos como el de Euler, el del Trapecio, el de Runge-Kutta y el de Heun.

Plasmaremos las soluciones en gráficos que nos permitirán dar una interpretación físico química a los resultados.

2 Reacción I

En primer lugar consideramos una reacción química irreversible en una solución bien mezclada, que se desarrolla para un volumen y una temperatura constantes, y en la que uno de los componentes (B, en este caso), hace de efecto catalítico.

[math] A + B →_{k1} 2B [/math]

Para desarrollar nuestro análisis nos apoyaremos en dos leyes físico químicas:

Ley de acción de masas: Enunciada por Guldberg y Waage en 1864, establece, entre otras cosas, que la velocidad de reacción es proporcional al producto de las concentraciones de los reactivos.

Principio de conservación de la masa: Establece que en una reacción química ordinaria la masa permanece constante, siendo la masa consumida de los reactivos igual a la masa obtenida de los productos. Fue elaborada independientemente por Míjail Lomonósov en 1745 y Antoine Lavoisier en 1785.

2.1 Interpretación del problema y deducción de las ecuaciones diferenciales

Para hallar las ecuaciones diferenciales que nos describan cómo se desarrolla la reacción en el tiempo nos basaremos en las dos leyes enunciadas anteriormente. En primer lugar, definimos nuestra variable independiente como el tiempo [math]t[/math] en segundos. Las dos variables dependientes serán las concentraciones de los compuestos A y B, que denominaremos [math]x(t)[/math] e [math]y(t)[/math], ambas expresadas en mol/l.

El principio de conservación de la masa expone que la suma de las concentraciones de los reactivos y los productos es constante. Esto equivale, derivando la expresión [math]x(t)+y(t)=cte[/math], a que desaparece tanto reactivo como producto aparece, o lo que es lo mismo, que la velocidad de aparición del producto es igual a la velocidad de desaparición (velocidad de aparición negativa) del reactivo:

[math]x'(t)=-y'(t)[/math]

La otra ecuación que necesitamos la obtenemos directamente de aplicar la Ley de acción de masas, que establece que la velocidad de la reacción (la velocidad a la que aparece el producto o desaparece el reactivo) es igual al producto de las concentraciones por una constante.

[math]y'(t)=k_1x(t)y(t)[/math]

2.2 Obtención del problema de valor inicial

Una vez demostrado el origen de las ecuaciones (1) y (2) en el apartado anterior, obtendremos el problema de valor inicial (PVI) mediante las operaciones detalladas a continuación:

[math]x'(t)+y'(t)=0\quad\quad(1)[/math]

[math]y'(t)=k_1x(t)y(t)\quad\quad(2)[/math]

Integrando (1) obtenemos (3):

[math]\int_\mathbb{D} (\frac{dx(t)}{dt} + \frac{dy(t)}{dt})\,dt =\int_\mathbb{D} 0\,dt \quad \Rightarrow \quad x(t)+y(t)=k_2\quad \Rightarrow \quad x(t)=k_2-y(t)\quad\quad(3)[/math]

Sustituyendo (3) en (2), obtenemos (4):

[math]y'(t)=k_1x(t)y(t)\quad \Rightarrow \quad y'(t)=k_1(k_2-y(t))y(t)\quad \Rightarrow \quad y'(t)=k_1k_2y(t)-k_1y^2(t)\quad\quad(4)[/math]

Conocemos las concentraciones iniciales de A y de B, por lo que podemos obtener el valor de la constante k₂:

[math][A]:\quad x(0)=1 mol/l;\quad\quad [B]:\quad y(0)=0,01 mol/l[/math]

[math]x(t)+y(t)=k_2\quad \Rightarrow \quad x(0)+y(0)=k_2\quad \Rightarrow \quad 1+0,01=k_2\quad \Rightarrow \quad k_2=1,01[/math]

Por último, conocemos también el valor de k₁, por lo que podremos plantear a continuación el correspondiente PVI:

[math]k_1=1[/math]

[math](PVI)\begin{cases}y'(t)=1,01y(t)-y^2(t) \\ \\ y(0) = 0,01 \end{cases}\quad (t,y)\in[0,10]\times\mathbb{R}[/math]

2.3 Unicidad de la solución

Intersección B con dominio

Aplicamos a continuación el Teorema de Existencia y Unicidad:

Sea [math]f(t,y)=y'(t) \Rightarrow f(t,y)=1.01y(t)-y^2(t)\quad[/math], y [math]\quad\frac{\partial }{\partial y} f(t,y) = 1.01-2y(t)[/math].

Y sea [math]B[(t_0,y_0),r],\quad r\gt0,\quad en (t_0,y_0)=(0,0.01)[/math].

[math]f(t,y)[/math] es continua en [math] \mathbb{B}\cap \mathbb{D}\quad \Rightarrow \quad \exists y(t)\quad /\quad y'(t)=1,01y(t)-y^2(t)[/math]

[math]\frac{\partial }{\partial y} f(t,y)[/math] es continua en [math] \mathbb{B}\cap \mathbb{D}\quad \Rightarrow \quad \exists ! y(t)\quad /\quad y'(t)=1,01y(t)-y^2(t)[/math]

Con lo que queda demostrado que el (PVI) tiene solución y además esta es única.

2.4 Resolución numérica del problema de valor inicial

2.4.1 Método de Euler

[math] PVI = \left \{ \begin{array}{rcl} y_0\\ y_{n+1} = y_n + h*f(t_n,y_n)\\ \end{array} \right . [/math]

Representación gráfica y(t). Método de Euler

%Método de Euler
clc

%Datos iniciales
t0=0;
tN=4;
y0=0.01;
h=0.1;
t0=0;
tN=10;

%Cálculo de subintervalos
N=round((tN-t0)/h);
t=t0:h:tN;

%Creamos el vector donde vamos a guardar las soluciones
y=zeros(1,N+1);
y(1)=y0;

%Euler
for i=1:N;
    y(i+1)=y(i)+h*(1.01*y(i)-(y(i)).^2);
end
%Tabla de Resultados
[t',y']
x=1.01-y;
hold on 
plot(t,y,'r')
plot(t,x,'b')
hold off
legend('Sustancia B: y(t)','Sustancia A: x(t)','Location','best');
xlabel('Tiempo (s)')
ylabel('Concentración (mol/l)')
grid on

2.4.2 Método del trapecio

[math] PVI = \left \{ \begin{array}{rcl} y_0\\ y_{n+1} = y_n + \frac{h}{2}*f(t_n,y_n)\\ \end{array} \right . [/math]

Representación gráfica y(t). Método del Trapecio

%Método del Trapecio
clc

%Datos iniciales
t0=0;
tN=4;
y0=0.01;
h=0.1;
t0=0;
tN=10;

%Cálculo de subintervalos
N=round((tN-t0)/h);
t=t0:h:tN;

%Creamos el vector donde vamos a guardar las soluciones
y=zeros(1,N+1);
y(1)=y0;

for i=1:N
y(i+1)=(((1-(h/2)*1.01)-sqrt(((1-(h/2)*1.01)^2)+4*(h/2)*((1+(h/2)*1.01)*y(i)-(h/2)*(y(i))^2)))/-h);
end
%Sacamos la tabla de resultados
[t',y']
x=1.01-y;
%Gráfico
hold on
plot(t,y,'r')
plot(t,x,'b')
legend('Sustancia B: y(t)','Sustancia A: x(t)','Location','best');
xlabel('Tiempo (s)')
ylabel('Concentración (mol/l)')
grid on

2.4.3 Método de Runge-Kutta

El método de resolución numérica de Runge-Kutta, involucra un gran número de parámetros. Para simplificar la implementación de este método en Matlab, utilizaremos la función inline, con la que declararemos previamente la función algebraica necesaria para el cálculo de los parámetros [math]k_i[/math], ahorrándonos tener que escribir innumerables veces la correspondiente expresión matemática.

A continuación, detallamos la función [math]f(y)[/math] de la que dependerán los valores de los parámetros [math]k_i[/math]:

[math]f(y)=y'(t)=1.01y(t)-y(t)^2 \quad \Rightarrow \quad \begin{cases} k_1 & = f(y_i) \\ k_2 & = f(y_i + {1 \over 2}k_1h) \\ k_3 & = f(y_i + \frac{1}{2}k_2h) \\ k_4 & = f(y_i + k_3h) \\ \end{cases}\quad [h=0.1] [/math]

Entonces el método RK4 para este problema está dado por la siguiente ecuación:

[math]y_{i+1} = y_{i} + {1 \over 6}h(k_1+2k_2+2k_3+k_4)[/math]

Mediante el siguiente programa en Matlab calculamos numéricamente la solución y la representamos en un gráfico:

Representación gráfica y(t). Método Runge-Kutta

%Método de Runge-Kutta
clc

%Función f(t)=y'(t)
fy=inline('1.01*y-y^2','y');

%Datos del problema
y0=0.01;
t0=0;
tN=10;
h=0.1;

%Variable independiente (t)
t=t0:h:tN;

%Variable dependiente (y)
y=zeros(1,length(t));
y(1)=y0;

for i=1:length(t)-1
    k1=fy(y(i));
    k2=fy(y(i)+k1/2*h);
    k3=fy(y(i)+k2/2*h);
    k4=fy(y(i)+k3*h);
    y(i+1)=y(i)+h/6*(k1+2*k2+2*k3+k4);
end

%Variable dependiente (x)
x=1.01-y;

%Representación gráfica
hold on
plot(t,x,'b')
plot(t,y,'r')
legend('Sustancia A: x(t)','Sustancia B: y(t)')
xlabel('Tiempo (s)')
ylabel('Concentración (mol/l)')
grid on

2.4.4 Comparación de métodos

Comparación de métodos

Comparación de métodos

Como podemos comprobar en las imágenes, el método de Euler es menos preciso que el del trapecio o el de Runge-Kutta.

2.5 Resolución en forma de sistema de ecuaciones

A continuación vamos a plantear el problema de valor inicial con un sistema de dos ecuaciones equivalente al problema resuelto anteriormente.

Estas ecuaciones serán:

[math](PVI) \begin{cases} y'(t)=k_1 x(t)y(t) \\ x'(t)=-k_1 x(t)y(t) \\ \\ x(0)=1 \\ y(0) = 0.01 \end{cases} \quad t\in[0,10][/math]

Resolveremos el problema por el método de Euler y por el método de Runge-Kutta, plasmando las soluciones en gráficos.

2.5.1 Método de Euler

Representación gráfica de las concentraciones de A y B. Resolución del sistema por Euler

%Resolución del sistema por el método de Euler

clear all

%Definimos el vector de la variable independiente, el número de
%subintervalos y la longitud de paso

t0 = 0;
tN = 10;
h = 0.1;
t = t0:h:tN
N = (tN-t0)/h;

%Definimos los vectores de las variables dependientes

x = zeros(1, N+1);
y = zeros(1, N+1);

x(1) = 1;
y(1) = 0.01;

for i = 1:N
    
    x(i+1)=x(i)+h*(-x(i).*y(i));
    y(i+1)=y(i)+h*(y(i).*x(i));
    
end

%Representamos gráficamente las soluciones

hold on

plot(t,x,'b');
plot(t,y,'r');
grid on

legend('Sustancia A: x(t)','Sustancia B: y(t)')
xlabel('Tiempo (s)')
ylabel('Concentración (mol/l)')


hold off

2.5.2 Método de Runge-Kutta

El método de resolución numérica de Runge-Kutta, involucra un gran número de parámetros. Para simplificar la implementación de este método en Matlab, el sistema se tratará como una matriz [math]z[/math]. Además, utilizaremos la función inline, con la que declararemos previamente la función algebraica necesaria para el cálculo de los parámetros [math]k_i[/math], ahorrándonos tener que escribir innumerables veces la correspondiente expresión matemática.

Primero, expresamos el sistema en forma matricial:

[math](PVI)\begin{cases} x'(t)=-k_1 x(t)y(t) \\ y'(t)=k_1 x(t)y(t) \\ \\ x(0)=1 \\ y(0) = 0.01 \end{cases}\quad t\in[0,10] \quad \Rightarrow \quad k_1=1\quad \Rightarrow \quad (PVI)\begin{cases}\vec{z}'(t)= \begin{bmatrix} x'(t) \\ y'(t) \\ \end{bmatrix} = \begin{bmatrix} -x(t)y(t) \\ x(t)y(t) \\ \end{bmatrix} \\ \\ \vec{z}(0)= \begin{bmatrix} x(0) \\ y(0) \\ \end{bmatrix} = \begin{bmatrix} 1 \\ 0.01 \\ \end{bmatrix} \end{cases}\quad t\in[0,10][/math]

A continuación, detallamos la función [math]\vec{F}[/math] de la que dependerán los valores de los parámetros [math]k_i[/math]:

[math]\vec{F}(x,y)=\vec{F}(\vec{z})= \begin{bmatrix} f(\vec{z}) \\ g(\vec{z}) \\ \end{bmatrix}= \begin{bmatrix} x'(t) \\ y'(t) \\ \end{bmatrix}= \vec{z}'(t)= \begin{bmatrix} -x(t)y(t) \\ x(t)y(t) \\ \end{bmatrix} \quad \Rightarrow \quad \begin{cases} [\vec{k_1}] & = F(\vec{z}_i) \\ [\vec{k_2}] & = F(\vec{z}_i + {1 \over 2}k_1h) \\ [\vec{k_3}] & = F(\vec{z}_i + \frac{1}{2}k_2h) \\ [\vec{k_4}] & = F(\vec{z}_i + k_3h) \\ \end{cases}\quad [h=0.1] [/math]

Entonces el método RK4 para este problema está dado por la siguiente ecuación:

[math]\vec{z}_{i+1} = \vec{z}_{i} + {1 \over 6}h(\vec{k_1}+2\vec{k_2}+2\vec{k_3}+\vec{k_4})[/math]

Mediante el siguiente programa en Matlab calculamos numéricamente la solución y la representamos en un gráfico:

Representación gráfica de las concentraciones de A y B. Resolución del sistema por Runge-Kutta

%Método de Runge-Kutta
clc

%Función F=[x'(t); y'(t)]
F=inline('[-z(1)*z(2);z(1)*z(2)]','z');

%Datos del problema 
x0=1;
y0=0.01;
t0=0;
tN=10;
h=0.1;

%Vector de variable independiente (t)
t=t0:h:tN;

%Matriz de variables dependientes z=[x;y]
z=zeros(2,length(t));
z(:,1)=[x0;y0];

for i = 1:length(t)-1
    k1=F(z(:,i));
    k2=F(z(:,i)+k1.*h/2);
    k3=F(z(:,i)+k2.*h/2);
    k4=F(z(:,i)+k3.*h);
    z(:,i+1)=z(:,i)+h/6.*(k1+2.*k2+2.*k3+k4);
end
 
%Representación gráfica
hold on
plot(t,z(1,:))
plot(t,z(2,:),'r')
legend('Sustancia A: x(t)','Sustancia B: y(t)')
xlabel('Tiempo (s)')
ylabel('Concentración (mol/l)')
grid on
%hold off

2.6 Interpretación de las gráficas

Podemos observar por los resultados obtenidos que la velocidad de aparición de B es igual a la velocidad de desaparición de A (como dicta el principio de conservación de la masa), que A acaba desapareciendo prácticamente por completo y que la concentración de B termina por estabilizarse alrededor de 1 mol/l.

También cabe destacar que cuando el tiempo transcurrido alcanza 4,5 segundos las concentraciones de ambos compuestos son las mismas.

3 Reacción II

En este apartado vamos a estudiar la reacción propuesta por Alfred J. Lotka en 1920, que es la siguiente:

[math] A + X →_{k1} 2X;[/math]  :[math] X + Y →_{k2} 2Y;[/math]  :[math]Y →_{k3} 2B [/math]

Como se puede observar, se trata de una reacción compleja integrada por 3 reacciones simples (las dos primeras autocatalíticas) en las que intervienen un total de 4 compuestos, A, B, X e Y. Siguiendo el desarrollo de la reacción, es lógico suponer que A se agotará al comienzo de la reacción, que los máximos de concentración de X e Y se darán en los estados intermedios y que al final el compuesto B será el único restante, siendo el producto final de la reacción.

3.1 Interpretación del problema y deducción de las ecuaciones diferenciales

Para averiguar cómo evolucionan las concentraciones de los cuatro compuestos durante el desarrollo de la reacción necesitaremos obtener cuatro ecuaciones diferenciales. Para ello nos volveremos a basar en los principios físicos que ya empleamos en la anterior reacción. Definimos A(t), x(t), y(t), y B(t) (A, x, y, B; para ahorrar notación) a las concentraciones de A, X, Y, y B, respectivamente.

En primer lugar aplicamos el Principio de Conservación de la masa, que establece que la suma de las concentraciones de los compuestos es constante (no desaparece masa en el proceso), lo que derivando nos indica que:

[math] A' + y' + y' + B' = 0 \quad\quad(I)[/math]

A continuación aplicaremos la Ley de Acción de Masas (velocidad de reacción proporcional al producto de la concentración de los reactivos).

Aplicando esta ley a la última reacción simple obtenemos que:

[math]B'=k_3y \quad\quad(II)[/math]

En el caso de las concentraciones de X e Y las expresiones resultantes son algo más complicadas ya que estos dos compuestos actúan como reactivos y productos en más de una reacción simple, por lo que la variación de sus concentraciones será la suma de las variaciones de las reacciones en las que aparezcan:

[math] x' = x'_{1} + x'_{2} → x' = k_{1}Ax - k{2}xy \quad\quad(III)[/math]

[math] y' = y'_{2} + y'_{3} → y' = k_{2}xy − k_3{y} \quad\quad(IV)[/math]

Por último sustituimos (II), (III) y (IV) en (I) obteniendo:

[math] A' + k_{1}Ax - k{2}xy + k_{2}xy − k_3{y}+ k_3y= 0 : [/math]

De donde obtenemos que:

[math] A' = -k_{1}Ax [/math]

Ahora que ya tenemos las variaciones de concentración de los compuestos en función de las concentraciones podemos aplicar métodos numéricos que nos lleven a una solución aproximada.

3.2 Resolución numérica del sistema de ecuaciones

3.2.1 Método de Euler

Representación gráfica de las concentraciones de A, B, X, e Y a lo largo del tiempo. Resoluciónd delsistema por el método de Euler con longitud de paso h=0.01

Representación gráfica de las concentraciones de A, B, X, e Y a lo largo del tiempo. Resoluciónd delsistema por el método de Euler con longitud de paso h=0.001

%Resolución del sistema de la reacción de Lotka por el método de Euler

clear all

%Definimos la variable independiente y establecemos los datos del problema
%de valor inicial. Para la longitud de paso h introducimos en primer lugar
%0.01 y luego 0.001

t0=0;
tN=200;
h= input('Introduzca una longitud de paso ');
t = t0:h:tN;
N=(tN-t0)/h;
k1 = 0.1;
k2 = 0.1;
k3 = 0.05;
A0 = 5;
B0 = 0;
X0 = 5*10^(-4);
Y0 = 10^(-5);


%Definimos las variables dependientes.

A = zeros(1, N-1);
B = A;
X = A;
Y = A;

A(1) = A0;
B(1) = B0;
X(1) = X0;
Y(1) = Y0;

for i = 1:N
    
    X(i+1) = X(i) + h*(k1*A(i)*X(i)-k2*Y(i)*X(i));
    Y(i+1)=Y(i)+h.*(k2.*X(i).*Y(i)-k3.*Y(i));
    B(i+1) = B(i) + h*(k3*Y(i));
    A(i+1) = A(i) + h*(-k1*A(i)*X(i));
   
end

%Representamos gráficamente la solución

hold on

plot (t,A,'b')
plot(t, X, 'r')
plot(t,Y, 'g')
plot(t, B, 'm')
legend('Sustancia A: A(t)','Sustancia X: x(t)','Sustancia Y: y(t)','Sustancia B: B(t)','Location','Best')
xlabel('Tiempo (s)')
ylabel('Concentración (mol/l)')
grid on

hold off

3.2.2 Método de Heun

Si bien el método de resolución numérica de Heun no involucra un número de parámetros tan notable como el método de Runge-Kutta, en este caso tenemos un sistema de 4 ecuaciones diferenciales, por lo que, en aras de simplificar, utilizaremos el mismo criterio que en el apartado 2.5.2.: El sistema se tratará como una matriz [math]z[/math], y para la obtención de los parámetros [math]k_i[/math] utilizaremos la función inline, con la que declararemos previamente las expresiones necesarias para el cálculo de los parámetros [math]k_i[/math].

Primero, expresamos el sistema en forma matricial:

[math](PVI)\begin{cases} A'(t)=-k_1 A(t)x(t) \\ x'(t)=k_1 A(t)x(t)-k_2 x(t)y(t) \\ y'(t)=k_2 x(t)y(t)-k_3 y(t) \\ B'(t)=k_3 x(t)y(t) \\ \\ A(0)=5 \\ x(0)=5·10^{-4} \\ y(0)=10^{-5} \\ B(0) = 0 \end{cases}\quad t\in[0,200] \quad \Rightarrow \quad \begin{cases} k_1=0.1 \\ k_2=0.1 \\ k_3=0.05 \\ \end{cases} \quad \Rightarrow \quad (PVI)\begin{cases}\vec{z}'(t)= \begin{bmatrix} A'(t) \\ x'(t) \\ y'(t) \\ B'(t) \\ \end{bmatrix} = \begin{bmatrix} -0.1 A(t)x(t) \\ 0.1 A(t)x(t)-0.1 x(t)y(t) \\ 0.1 x(t)y(t)-0.05 y(t) \\ 0.05 y(t)x(t)\\ \end{bmatrix} \\ \\ \vec{z}(0)= \begin{bmatrix} A(0) \\ x(0) \\ y(0) \\ B(0) \\ \end{bmatrix} = \begin{bmatrix} 5 \\ 5·10^{-4} \\ 10^{-5} \\ 0 \end{bmatrix} \end{cases}\quad t\in[0,200][/math]

A continuación, detallamos la función [math]\vec{F}[/math] de la que dependerán los valores de los parámetros [math]k_i[/math]:

[math]\vec{F}(A,x,y,B)=\vec{F}(\vec{z})= \begin{bmatrix} f(\vec{z}) \\ g(\vec{z}) \\ h(\vec{z}) \\ i(\vec{z}) \\ \end{bmatrix}= \begin{bmatrix} A'(t) \\ x'(t) \\ y'(t) \\ B'(t) \\ \end{bmatrix}= \vec{z}'(t)= \begin{bmatrix} -0.1 A(t)x(t) \\ 0.1 A(t)x(t)-0.1 x(t)y(t) \\ 0.1 x(t)y(t)-0.05 y(t) \\ 0.05 y(t)x(t)\\ \end{bmatrix} \quad \Rightarrow \quad \begin{cases} [\vec{k_1}] & = F(\vec{z}_i) \\ [\vec{k_2}] & = F(\vec{z}_i + {1 \over 2}k_1h) \\ \end{cases}\quad [h=0.01] [/math]

Entonces el método de Heun para este problema está dado por la siguiente ecuación:

[math]\vec{z}_{i+1} = \vec{z}_{i} + {1 \over 2}h(\vec{k_1}+\vec{k_2})[/math]

Mediante el siguiente programa en Matlab calculamos numéricamente la solución y la representamos en un gráfico:

Representación gráfica de las concentraciones de A,X,Y y B. Resolución del sistema por Heun

%Método de Heun
clc

%Función F=[A'(t);x'(t); y'(t);B'(t)]
F=inline('0.1.*[-z(1)*z(2);z(1)*z(2)-z(2)*z(3);z(2)*z(3)-z(3)/2;z(3)/2]','z');

%Datos del problema
A0=5;
x0=5*10^-4;
y0=10^(-5);
B0=0;
t0=0;
tN=200;
h=0.1;

%Vector de variable independiente (t)
t=t0:h:tN;

%Matriz de variables dependientes z=[A;x;y;B]
z=zeros(4,length(t));
z(:,1)=[A0;x0;y0;B0];

for i = 1:length(t)-1
    k1=F(z(:,i));
    k2=F(z(:,i)+k1.*h/2);
    z(:,i+1)=z(:,i)+h/2.*(k1+k2);
end
 
%Representación gráfica
hold on
plot(t,z(1,:),'b')
plot(t,z(2,:),'r')
plot(t,z(3,:),'g')
plot(t,z(4,:),'m')
legend('Sustancia A: A(t)','Sustancia X: x(t)','Sustancia Y: y(t)','Sustancia B: B(t)','Location','Best')
xlabel('Tiempo (s)')
ylabel('Concentración (mol/l)')
grid on
hold off

3.3 Interpretación de las gráficas