Figure 5.7 (page 210):

Full model solution for k_1=1, k_{-1}=0.5, k_2=k_{-2}=10.

Code for Figure 5.7

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 k1 k_1 k2 k_2 nscheme
k1      = 1;
k_1     = 1/2;
k2      = 10;
k_2     = 10;
c_A0    = 1;
c_B0    = 0.4;
c_C0    = 0;
%
% nscheme: 0 full model, 1 1st is rls, 2 2nd is rls; 3 qssa on B 
%
nscheme = 0;


if (nscheme == 0)
  %
  % IC for full model
  %
  x0=[c_A0; c_B0; c_C0];
elseif (nscheme == 1)
  %
  % IC for rx 1 is rls
  %
  K2    = k2/k_2;
  x0=[c_A0; 1/(1+K2)*(c_B0+c_C0); K2/(1+K2)*(c_B0+c_C0)];
elseif (nscheme == 2)
  %
  % IC for rx 2 is rls
  %
  K1    = k1/k_1;
  x0=[1/(1+K1)*(c_A0+c_B0); K1/(1+K1)*(c_A0+c_B0); c_C0];
elseif (nscheme == 3)
  %
  % IC for QSSA model
  %
  x0=[c_A0; k1/(k_1+k2)*c_A0 + k_2/(k_1+k2)*c_C0; c_C0];
else
  error ('nscheme out of range');
end

tfinal = 5;
ntimes = 200;
tout   = linspace(0, tfinal, ntimes)';
opts = odeset ('AbsTol', sqrt (eps), 'RelTol', sqrt (eps));
[tsolver, x] = ode15s (@rhs, tout, x0, opts);
%
% compute RA RB RC production rates
%
RA = -k1*x(:,1)+k_1*x(:,2);
RC = k2*x(:,2)-k_2*x(:,3);
RB = -RA -RC;

table = [tout x RA RB RC];


plot (table(:,1),table(:,2:4));
title ('Figure 5.7')

rhs.m

function retval = rhs(t,x)
global k1 k_1 k2 k_2 nscheme
c_A   = x(1);
c_B   = x(2);
c_C   = x(3);
r1  = k1*c_A - k_1*c_B;
r2  = k2*c_B - k_2*c_C;
retval = zeros (3, 1);
if (nscheme == 0)
   %
   % full model
   %
   retval(1) = -r1; 
   retval(2) = r1 - r2;
   retval(3) = r2;
elseif (nscheme == 1)
   %
   % reaction 1 is rls, reaction 2 at equilibrium
   %
   K2    = k2/k_2;
   retval(1) = -r1; 
   retval(2) =  1/(1+K2)*r1;
   retval(3) =  K2/(1+K2)*r1;
elseif (nscheme == 2)
  %
  % reaction 2 is rls, reaction 1 at equilibrium
  %
  K1    = k1/k_1;
  retval(1) = -1/(1+K1) * r2; 
  retval(2) = -K1/(1+K1) * r2;
  retval(3) = r2;
elseif (nscheme == 3)
  %
  % QSSA on B, redefine the c_B concentration under QSSA
  %
  c_B= k1/(k_1+k2)*c_A + k_2/(k_1+k2)*c_C;
  r1  = k1*c_A - k_1*c_B;
  r2  = k2*c_B - k_2*c_C;
  retval(1) = -r1; 
  %retval(2) = r1 - r2; don't use this one, it is always zero and
  %therefore incorrect
  retval(2) = k1*k2*(k_2-k1)/(k_1+k2)^2 * c_A +  ...
              k_1*k_2*(k1-k_2)/(k_1+k2)^2*c_C;
  retval(3) = r2;
else
  error ('nscheme out of range');
end