function BatchReactor
% 間歇反應器中的連串-平行復雜反應系統
%
%   Author: Huang Huajiang
%   Copyright 2003 UNILAB Research Center, 
%   East China University of Science and Technology, Shanghai, PRC
%   $Revision: 1.0 $  $Date: 2002/07/16 $

clear all
clc

T = 224.6 + 273.15; % Reactor temperature, Kelvin
R = 8.31434;	    % Gas constant, kJ/kmol K

% Arrhenius constant, 1/s
k0 = [5.78052E+10  3.92317E+12  1.64254E+4  6.264E+8];

% Activation energy, kJ/kmol
Ea = [124670  150386  77954  111528];

% 初始濃度C0(i), kmol/m^3
C0 = [1  0  0  0  0];
tspan = [0 1e4];
[t,C] = ode45(@MassEquations, tspan, C0,[],k0,Ea,R,T)

% 繪圖
plot(t,C(:,1),'r-',t,C(:,2),'k:',t,C(:,3),'b-.',t,C(:,4),'k--');
xlabel('Time (s)');
ylabel('Concentration (kmol/m^3)');
legend('A','B','C','D')
CBmax = max(C(:,2));    % CBmax: the maximum concentration of B, kmol/m^3
yBmax = CBmax/C0(1)     % yBmax: the maximum yield of B
index = find(C(:,2)==CBmax);
t_opt = t(index)        % t_opt: the optimum batch time, s

% ------------------------------------------------------------------
function dCdt = MassEquations(t,C,k0,Ea,R,T)
% Reaction rate constants, 1/s
k = k0.*exp(-Ea/(R*T));
k(5) = 2.16667E-04;  

% Reaction rates, kmoles/m3 s
rA = -(k(1)+k(2))*C(1);
rB = k(1)*C(1)-k(3)*C(2);
rC = k(2)*C(1)-k(4)*C(3);
rD = k(3)*C(2)-k(5)*C(4);
rE = k(4)*C(3)+k(5)*C(4);

% Mass balances
dCdt = [rA; rB; rC; rD; rE];