function d_depnt__d_indep = pfr(indep, depnt)
% Dynamic balance for the reactor
%
% indep: the independent ODE variable, such as time or length
% depnt: a vector of dependent variables
%
% X = depnt(1) = the conversion
% y = depnt(2) = the pressure ratio = P/P_0 = y
%
% Returns d(depnt)/d(indep) = a vector of ODEs
% Assign some variables for convenience of notation
FA = depnt(1);
FB = depnt(2);
FC = depnt(3);
FD = depnt(4);
FE = depnt(5);
FG = depnt(6);
FW = depnt(7);
y = depnt(8);
% Constant(s)
k1 = 0.014; % L^{0.5} / mol^{0.5} / s
k2 = 0.007; % L/(mol.s)
k3 = 0.14; % 1/s
k4 = 0.45; % L/(mol.s)
alpha = 0.002; % 1/L
CT0 = 1.0; % mol/L
FA0 = 10; % mol/s
FB0 = 5.0; % mol/s
FT0 = FA0 + FB0;
FT = FA + FB + FC + FD + FE + FW + FG;
CA = CT0 * FA/FT * y;
CB = CT0 * FB/FT * y;
CC = CT0 * FC/FT * y;
CD = CT0 * FD/FT * y;
CE = CT0 * FE/FT * y;
CG = CT0 * FG/FT * y;
CW = CT0 * FW/FT * y;
% Reaction 1: A + 0.5B -> C
r1A = -k1*(CA)*(CB)^(0.5);
r1B = 0.5*r1A;
r1C = -r1A;
%# Reaction 2: 2A -> D
r2A = -k2*(CA)^2;
r2D = -r2A/2;
% Reaction 3: C -> E + W
r3C = -k3*(CC);
r3E = -r3C;
r3W = -r3C;
% Reaction 4: D + W -> G + C
r4D = -k4*(CD)*(CW);
r4W = r4D;
r4G = -r4D;
r4C = -r4D;
% Output from this ODE function must be a COLUMN vector, with n rows
n = numel(depnt);
d_depnt__d_indep = zeros(n,1);
d_depnt__d_indep(1) = r1A + r2A;
d_depnt__d_indep(2) = r1B;
d_depnt__d_indep(3) = r1C + r3C + r4C;
d_depnt__d_indep(4) = r2D + r4D;
d_depnt__d_indep(5) = r3E;
d_depnt__d_indep(6) = r3W + r4W;
d_depnt__d_indep(7) = r4G;
d_depnt__d_indep(8) = -alpha/(2*y) * (FT / FT0);
% Integrate the ODE
% -----------------
% The independent variable always requires an initial and final value:
indep_start = 0.0; % kg
indep_final = 500.0; % kg
% Set initial condition(s): for integrating variables (dependent variables)
FA_depnt_zero = 10.0; % i.e. FA(W=0) = 10.0
FB_depnt_zero = 5.0; % i.e. FB(W=0) = 10.0
FC_depnt_zero = 0.0; % i.e. FC(W=0) = 10.0
FD_depnt_zero = 0.0; % etc
FE_depnt_zero = 0.0;
FG_depnt_zero = 0.0;
FW_depnt_zero = 0.0;
y_depnt_zero = 1.0; % i.e. y(W=0) = 1.0
IC = [FA_depnt_zero, FB_depnt_zero, FC_depnt_zero, FD_depnt_zero, ...
FE_depnt_zero FG_depnt_zero, FW_depnt_zero, y_depnt_zero];
% Integrate the ODE(s):
[indep, depnt] = ode45(@pfr, [indep_start, indep_final], IC);
[indep_p, depnt_p] = ode45(@pfr_pressure, [indep_start, indep_final], IC);
% Calculate Yields and Selectivities
FA = depnt(:,1);
FC = depnt(:,3);
FD = depnt(:,4);
FE = depnt(:,5);
Yield_C = FC ./ (FA_depnt_zero - FA);
S_CE = FC./FE;
S_CD = FC./FD;
% Plot the results:
clf;
plot(indep, depnt(:,1), indep, depnt(:,2), indep, depnt(:,3), ...
indep, depnt(:,4), indep, depnt(:,5), indep, depnt(:,6), indep, depnt(:,7))
grid('on')
hold('on')
plot(indep, depnt(:,2), 'g')
xlabel('Catalyst weight, W [kg]')
ylabel('X and y')
legend('FA', 'FB', 'FC', 'FD', 'FE', 'FG', 'FW')
title('No pressure drops')
figure;
plot(indep, depnt(:,3), indep_p, depnt_p(:,3))
grid('on')
hold('on')
xlabel('Catalyst weight, W [kg]')
ylabel('X and y')
legend('FC (no pressure drop)', 'FC (with pressure drop)', 'Location', 'best')
figure
plot(indep_p, depnt_p(:,8))
grid('on')
xlabel('Catalyst weight, W [kg]')
ylabel('Pressure drop')