Figure 7.20 (page 392):

Dimensionless concentration versus radius for the nonisothermal spherical pellet: lower (A), unstable middle (B), and upper (C) steady states.

Code for Figure 7.20

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 Gamma beta Phi rout
%%
%% compute the pellet temperature and concentration profiles 
%% for the multiple steady state Weisz Hicks problem  
%%
%% 7/18/01
%%
%% jbr
%%
Gamma = 30;
beta  = 0.4;
Phiscale  = 0.01;

tol = 1e-6;
[value, nfun] = quadl (@integrand, 0, 1, tol);
intercept = sqrt(2*value);
Phi = intercept*Phiscale;

rout = linspace(0,3,100)';
%c0=1e-11;
c0vec = [1e-11; 0.5; 0.95];
results = [rout];
for i = 1: length(c0vec)
  c0 = c0vec(i);
  x0=[c0; 0; 1; 0];
  ode_opts = odeset ();
  rel = ode_opts.RelTol;
  if (isempty (rel))
    rel = sqrt (eps);
  end
  ode_opts = odeset('AbsTol',rel *c0,'RelTol',sqrt (eps));
  [rsolver,xout] = ode15s(@pelletode, rout, x0, ode_opts);
  %%
  %% checked sensitivies by finite difference, done 7/18/01
  %%
  ceinit = xout(length(rout),1);
  c0init = c0;
  y0     = [ ceinit; c0init];
  ydot0  = [-1; -1/xout(length(rout),3)];
  tsteps = linspace(0, ceinit-1 ,10)';
  dae_opts = odeset ('AbsTol', rel*c0init, 'RelTol', sqrt (eps));
  [tsolver,y]  = ode15i (@continode, tsteps, y0, ydot0, dae_opts);
  %%
  %% solve the pellet problem from correct c0
  %%
  c0 = y(length(tsteps),2);
  x0=[c0; 0; 1; 0];
  ode_opts = odeset ('AbsTol',rel*c0,'Reltol',sqrt (eps));
  [rsolver,xout] = ode15s(@pelletode, rout, x0, ode_opts);
  conc = xout(:,1);
  temp = beta*(1-conc);
  eta(i) = xout(length(rout),2)/Phi^2;
  results = [results conc temp];
end
subplot (2, 1, 1);
plot (results(:,1), results(:,2:2:6));
title ('Figure 7.20')

subplot (2, 1, 2);
plot (results(:,1), results(:,3:2:7));
title ('Figure 7.20')

integrand.m

function retval = integrand (c)
  global beta Gamma
  retval = c *exp(Gamma*beta*(1-c)/(1+beta*(1-c)));

pelletode.m

function xdot = pelletode(r,x)
  global Gamma beta Phi
  ca    = x(1);
  dcadr = x(2);
  s1    = x(3);
  s2    = x(4);
  xdot = zeros (4, 1);
  xdot(1) = dcadr;
  if (r == 0)
    xdot(2) = 1/3* Phi^2*ca*exp(Gamma*beta*(1-ca)/(1+beta*(1-ca)));
    xdot(4) = 0;
  else
    tmp = Phi^2*ca*exp(Gamma*beta*(1-ca)/(1+beta*(1-ca)));
    xdot(2) = - 2/r*dcadr + tmp;
    xdot(4) = tmp*(1/ca - Gamma*beta/(1+beta*(1-ca))^2)*s1 - 2/r*s2;
  end
  xdot(3) = s2;

continode.m

function res = continode(t, y, ydot)
  global rout
  ce = y(1);
  c0 = y(2);
  %%
  %% solve pellet problem from c0
  %%
  x0=[c0; 0; 1; 0];
  ode_opts = odeset('AbsTol',1e-5*c0, 'RelTol', sqrt (eps));
  [rsolver,xout] = ode15s(@pelletode, rout, x0, ode_opts);
  %%
  res = zeros (2, 1);
  res(1) = ydot(1) + 1;
  res(2) = ydot(2) + 1/xout(length(rout),3);