Figure A.3 (page 650):

Solution to first-order differential equation dc_A/dt=-kc_A, and sensitivities S_1=d c_A/d k and S_2=d c_A/d c_{A0}.

Code for Figure A.3

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 ca0 k
%%
%% Illustrate sensitivity calculation in simple first-order reaction
%%
%% dca/dt = -k ca
%% ca(0)  = ca0
%% 
%% theta = [k; ca0]
%% 
%% compute solution to model
%% and sensitivities 
%%
%%  ca(t)
%%  S1 = d ca / d k
%%  S2 = d ca / d ca0
%% 
%%  ca(t) = ca0 exp(-kt)
%%  S1 = -t ca0 exp(-kt)
%%  S2 = exp(-kt)
%%
%% jbr,  10/5/01
%%
%% converted from ddasac to cvodes, 3/31/07, jbr
%%

if (~ exist ('CVode'))
  if (exist ('startup_STB'))
    %%startup_STB;
    startup_STB ("/usr/share/octave/packages/3.2/sundialsTB");
  else
    error ('%s: requires the Sundials Toolbox, please see the README file for installation instructions');
  end
end

ca0 = 2;
k   = 1;
tfinal = 5;
nts    = 100;
tout = linspace(0,tfinal,nts)';

%%
%% initialize state and parameters
%%
x0 = ca0;
data.theta = [k, ca0];
%%
%% initialize  sensitivites, sx
sx0 = [0, 1];
%% sensitivities for which parameters
sensflag = [1;2];

[x, sx, status] = cvode_sens(@firstorder, [], x0, sensflag, data, sx0, tout);

ca = x;
dcadk   = sx(1,1,:);
dcadca0 = sx(1,2,:);
table = [tout, x(:), dcadk(:), dcadca0(:)];
plot (table(:,1), table(:,2:4));
title ('Figure A.3')

firstorder.m

function [xdot, flag, new_data] = firstorder(t, x, data)
  global ca0 k
  ca      = x(1);
  k       = data.theta(1);
  xdot = -k*ca;
  flag = 0;
  new_data = [];

cvode_sens.m