function IsothermCSTRs    
% 三個串聯的CSTR等溫反應器
%
%   Author: HUANG Huajiang
%   Copyright 2003 UNILAB Research Center, 
%   East China University of Science and Technology, Shanghai, PRC
%   $Revision: 1.0 $  $Date: 2002/07/20 $

clear all
clc

CA0 = 1.8;      % kmol/m^3
CA10 = 0.4;     % kmol/m^3
CA20 = 0.2;     % kmol/m^3
CA30 = 0.1;     % kmol/m^3
k = 0.5;        % 1/min
tau = 2;
stoptime = 2.9; % min
[t,y] = ode45(@Equations,[0 stoptime],[CA10 CA20 CA30],[],k,CA0,tau);
disp(' Results:')
disp('    t         CA1       CA2       CA3')
disp([t,y])
plot(t,y(:,1),'k--',t,y(:,2),'b:',t,y(:,3),'r-')
legend('CA_1','CA_2','CA_3')
xlabel('Time (min)')
ylabel('Concentration')

% ------------------------------------------------------------------
function dydt = Equations(t,y,k,CA0,tau)
CA1 = y(1);
CA2 = y(2);
CA3 = y(3);
dCA1dt = (CA0-CA1)/tau - k*CA1;
dCA2dt = (CA1-CA2)/tau - k*CA2;
dCA3dt = (CA2-CA3)/tau - k*CA3;
dydt = [dCA1dt; dCA2dt; dCA3dt];