function conc=conc(C)
    A=1; %m2
    F1=3; %m3/s
    F2=0.5*C(1)/A; %m3/s
    CA1=1.5; %mol/m3
    CB1=1; %mol/m3
    k=0.5;
    conc(1)=F1-F2;
    conc(2)=(F1*CA1-F2*C(2)-C(2)*(F1-F2)-k*C(2)*C(3)*C(1))/C(1);
    conc(3)=(F1*CB1-F2*C(3)-C(3)*(F1-F2)-k*C(2)*C(3)*C(1))/C(1);
    conc(4)=(-F2*C(4)+k*C(2)*C(3)*C(1)-C(4)*(F1-F2))/C(1);
    conc(5)=(-F2*C(5)+k*C(2)*C(3)*C(1)-C(5)*(F1-F2))/C(1);
endfunction

t=0:1:30;
V0=1.5; %m3
CA0=2.5; %mol/m3
CB0=0; %mol/m3
CC0=0; %mol/m3
CD0=0; %mol/m3

res=lsode("conc",[V0,CA0,CB0,CC0,CD0],t);

figure(1)
plot(t,res(:,2),"linewidth",1.5)
hold on
plot(t,res(:,3),"linewidth",1.5)
plot(t,res(:,4),"linewidth",1.5)
plot(t,res(:,5),"linewidth",1.5)
xlabel("t, s","fontsize",12)
ylabel("C_A, C_B, C_C, C_D, mol/m^3","fontsize",12)
leyenda=legend("C_A","C_B","C_C","C_D");
set(leyenda,"fontsize",12)

figure(2)
plot(t,res(:,1),"linewidth",1.5)
xlabel("t, s","fontsize",12)
ylabel("V, m^3","fontsize",12)
leyenda2=legend("V");
set(leyenda2,"fontsize",12)