clc
% this problem statement is taken from:
%
% "Dynamic Models in Chemistry: A Workshop of Computer Simulations
%    Using Electronic Spreadsheets"
% by Daniel E. Atkinson, Douglas C. Brower, Ronald W. McClard
% Edited: David S. Barkley
% Published: N. Simonson & Company, Marina del Rey (310) 301-2847
%
% Section 7.6 Reaction Rates, p. 3 (Model 7.13)
%
%            k/1
%           ---->
%  CO + H O        CO  + H
%        2  <----    2    2
%           k/-1
%
% d[CO]   d[H2O]
% ----- = ------ = -k/1 [CO][H2O]
%  dt       dt
%
% d[CO2]   d[H2]
% ------ = ----- = -k/-1 [CO2][H2]
%   dt      dt
%
% initial conditions CO=1 H2O=1 CO2=0 H2=0
%
% y(1)=CO y(2)=H2O y(3)=CO2 y(4)=H2
%
%
% share the reaction rate between 'driver' and derivative routine
k1 = 5; km1 = 1;
global k1 km1
%
%
disp('the m file with derivatives ');
type chemreac

clc

t0 = 0;
tfinal=input('enter final value of time (e.g. 1.2) ');
y0 = [1; 1; 0; 0];           % Define initial conditions.

tol = 1.e-3;		% Accuracy
trace = 1;
[t,y] = ode23('chemreac',t0,tfinal,y0,tol,trace);
plot(t,y), title('Gas-Phase Reaction CO+H2O <=> CO2 + H2'), pause