clear all; close all; clc;
% REACTOR RCAI NO ISOTÉRMICO

% parámetros
tau	= 1; % hr , V/v
ko = 9703*3600; % hr-1
mDH	= 5960; % kcal/kmol
Ea	= 11843; % kcal/kmol
rocp	= 500; % kcal/(m3ºC)
Tin	= 25+273; % K
CAin	= 10; % kmol/m3
UA_V	= 150; % kcal/(m3ºC hr)
Tc	= 25+273; % K
R = 1.987; % kcal/(kmol.K)

function f = rcai(y,t,par)
  CA = y(1); T = y(2);
  tau=par(1);ko=par(2); mDH=par(3); Ea=par(4); rocp=par(5);
  Tin=par(6); CAin=par(7); UA_V=par(8); Tc=par(9); R=par(10);
  f(1) = (CAin - CA)/tau - ko*exp(-Ea/R./T).*CA;
  f(2) = (Tin - T)/tau +(mDH)/rocp*ko*exp(-Ea/R./T).*CA-UA_V/rocp*(T-Tc);
endfunction

par= [tau,ko, mDH,Ea, rocp, Tin, CAin, UA_V, Tc,R];
CA0 = input('punto inicial, CA (kmol/m3) =  ');
T0 = input('punto inicial, T (K) =  ');
y0 = [CA0; T0];
t0 = 0;
t = 0:.1:10;

y = lsode(@(x,t) rcai(x,t,par),y0,t);

figure()
subplot(211)
plot(t,y(:,1))
xlabel('t (hr)')
ylabel('CA (kmol/m3)')
subplot(212)
plot(t,y(:,2))
xlabel('t (hr)')
ylabel('T (K)')