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

figure()
plot([2.359 8.5636],[368.1 311.2],'bx','markersize',10) % puntos estables
hold on
plot([5.518],[339.1],'ro','markersize',10)  % punto inestable
axis([0,15,250,400])
xlabel('CA (kmol/m3)')
ylabel('T (K)')

par= [tau,ko, mDH,Ea, rocp, Tin, CAin, UA_V, Tc,R];
t = 0:.1:20;

for i = 1:20
  [x10,x20] = ginput(1);
  y0=[x10,x20] ;
  y = lsode(@(x,t) rcai(x,t,par),y0,t);
  plot(y(:,1),y(:,2),'k')
end