clc; clear all; close all;
pkg load control;
% RCAI con reacción de segundo orden A --> B
k = 1; % L/mol/min
V = 1; % L
v = 1;  % L/min

Cin = 1; % mol/L

function dCdt = Cprima(t,C,par)
  k=par(1); v=par(2); V=par(3); Cin=par(4);
  dCdt = v/V*Cin - v/V*C - k*C^2;
endfunction

% buscamos primero el punto de estado estacionario
Co = 0; % mol/L
Cest = fsolve(@(x) Cprima(0,x,[k,v,V,Cin]), Co);

C0 = 0; % mol/L
t = 0:0.1:5; % min
C = lsode(@(C,t) Cprima(t,C,[k,v,V, Cin]), C0,t);

plot(t,C)
xlabel('t (min)')
ylabel('C (mol/L)')
pause
% otra forma de resolverlo es usando la función de transferencia
% pero en este caso introducimos una aproximación en la linealización

pause
tau = 1/(v/V+2*k*Cest); k = 1/(1+2*V/v*k*Cest);
s =tf('s');
g = k/(tau*s+1);

Du = 1; % mol/L, altura del escalón
y=Du*step(g,t);
C = y + C0;

hold on % mantiene la curva en la ventana de gráfico
plot(t,C, 'r')
legend('No lineal', 'Linealizada')
xlabel('t (min)')
ylabel('C (mol/L)')