Figure 7.26 (page 404):

Molar flow of A versus reactor volume for second-order, isothermal reaction in a fixed-bed reactor.

Code for Figure 7.26

Text of the GNU GPL.

main.m

%% Copyright (C) 2001, James B. Rawlings and John G. Ekerdt
%%
%% This program is free software; you can redistribute it and/or
%% modify it under the terms of the GNU General Public License as
%% published by the Free Software Foundation; either version 2, or (at
%% your option) any later version.
%%
%% This program is distributed in the hope that it will be useful, but
%% WITHOUT ANY WARRANTY; without even the implied warranty of
%% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
%% General Public License for more details.
%%
%% You should have received a copy of the GNU General Public License
%% along with this program; see the file COPYING.  If not, write to
%% the Free Software Foundation, 59 Temple Place - Suite 330, Boston,
%% MA 02111-1307, USA.

global Naf P Rg T rhop rhob Da k Rp n xa

%%
%% units: g, mol, sec, cm, K (but atm pressure)
%%

k = 2.25e5;
Naf = 10;
T = 550;
rhop = 0.68;
rhob = 0.60;
Da = 0.008;
Rp = 0.45;
Rg =  82.06;
P  =  4;
n  = 2;
xa = 0.75;

Vanal = -2/(-n+1)*(2^((n+1)/2))*( (1-xa)^((-n+1)/2) - 1 )* ...
         Naf * Rp/3 * sqrt( (n+1)/2/Da) /                  ...
         ( rhob/rhop*sqrt(k)*((P/(Rg*T))^((n+1)/2)) );

Vanal2 = 4*( (1-xa)^(-1/2) - 1 )* ...
         Naf * Rp/3 * sqrt(3/Da) /  ...
         ( rhob/rhop*sqrt(k)*((P/(Rg*T))^(3/2)) );





npts   = 100;
vfinal = 5e5;
vtotal = vfinal*linspace(0,1,npts)';
v0     = vtotal(1);
vsteps = vtotal;
x0=Naf;
opts = odeset ('Events', @stop, 'AbsTol', sqrt (eps), 'RelTol', sqrt (eps));
[vout,x] = ode15s (@pbr, vsteps, x0, opts);
if ( length(vout) == npts )
  fprintf ('hey, did not reach final conversion, increase rstop\n');
end
xplot1 = x;
vplot = vout/1000.;
VR = vplot(length(vplot));
opts = odeset ('Events', @stop, 'AbsTol', sqrt (eps), 'RelTol', sqrt (eps));
[voutasy,xasy] = ode15s (@pbrasy, vsteps, x0, opts);
if ( length(voutasy) == npts )
  fprintf ('hey, did not reach final conversion, increase rstop\n');
end
xasyplot = xasy;
vasyplot = voutasy/1000.;
VRasy = vasyplot(length(vasyplot));
%ca = P/(Rg*T) * xplot1/(2*Naf);
%Phi = Rp/3*sqrt((n+1)/2*k*ca/Da);
%eta = 1.0 ./ Phi .* ( 1.0 ./tanh(3*Phi) - 1.0 ./ (3*Phi) );
table = [vplot xplot1];
tableasy =  [vasyplot xasyplot];
Naout = (1-xa)*Naf;
Natop = (1-xa+0.10)*Naf;
dashedlines = ...
[0,      Naout, VRasy, 0,     VR, 0,     VR, 2, VRasy, 2;
 1.1*VR, Naout, VRasy, Natop, VR, Natop, VR, 3, VRasy, 3];
plot (table(:,1), table(:,2), ...
       tableasy(:,1), tableasy(:,2), ...
       dashedlines(:,1), dashedlines(:,2), ...
       dashedlines(:,3), dashedlines(:,4), ...
       dashedlines(:,5), dashedlines(:,6));
title ('Figure 7.26')

stop.m

function [retval, isterminal, direction] = stop(t,x)
  global Naf xa
  Na = x(1);
  retval = Na - (1-xa)*Naf;

  isterminal = 1;
  direction = 0;

pbr.m

function xdot = pbr (t,x)
  global P Rg T Naf rhop rhob Da k Rp n
  Na = x(1);
  ca = P/(Rg*T) * Na/(2*Naf);
  Phi = Rp/3*sqrt((n+1)/2*k*ca/Da);
  %%
  %% try either the approximate first-order solution or
  %% the large Phi asymptotic solution, eta=1/Phi
  %%
  eta = 1./Phi*( 1./tanh(3*Phi) - 1/(3*Phi) );
%  eta = 1./Phi;
  xdot = -rhob/rhop*eta*k*ca*ca;

pbrasy.m

function xdot = pbrasy (t,x)
  global P Rg T Naf rhop rhob Da k Rp n
  Na = x(1);
  ca = P/(Rg*T) * Na/(2*Naf);
  Phi = Rp/3*sqrt((n+1)/2*k*ca/Da);
  %%
  %% try either the approximate first-order solution or
  %% the large Phi asymptotic solution, eta=1/Phi
  %%
%  eta = 1./Phi*( 1./tanh(3*Phi) - 1/(3*Phi) );
  eta = 1./Phi;
  xdot = -rhob/rhop*eta*k*ca*ca;