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

% vamos a suponer que el reactor está funcionando en estado estacionario
% con una entrada de 1 mol/L y súbitamente (a t = 1 min) esta cambia a 1.5 mol/L 
Cin0 = 1; % mol/L

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

% buscamos primero el punto de estado estacionario
Co = 0.2; % mol/L, un valor cualquiera para empezar a iterar
Cs0 = fsolve(@(x) Cprima(0,x,[k,v,V,Cin0]), Co)

% a t = 1 cambio Cin de 1.0 a 1.5 mol/L
t = 1:0.1:5; % min
Cin1 = 1.5; % mol/L
C = lsode(@(C,t) Cprima(t,C,[k,v,V,Cin1]), Cs0,t);

plot([0 t],[Cs0 C'])
xlabel('t (min)');ylabel('C (mol/L)');title('usando lsode')

% nuevo punto de estado estacionario:
Co = Cs0; % mol/L, un valor cualquiera para empezar a iterar
Cs1 = fsolve(@(x) Cprima(0,x,[k,v,V,Cin1]), Co)

pause

% otra forma de resolverlo es usando la función de transferencia
tau = 1/(v/V+k); k = 1/(1+V/v*k);
s = tf('s');
g = k/(tau*s+1);

Du = 0.5; % mol/L, altura del escalón (de 1.0 a 1.5 mol/L)
t1 = 0:0.1:4; % el comando step está pensado para considerar el tiempo a partir de cero
y = Du*step(g,t1); % step es un escalón unitario; variables desviación
C1 = y + Cs0;
figure(2)
plot([0 t],[Cs0 C1'])
xlabel('t (min)');ylabel('C (mol/L)');title('usando g(s) y step')