clear
reset(0)
echo on
format compact
%In this example we consider the performance
%of a batch reactor.  We have the sequential
%second order irreversible interactions:
%
%  A + B -> C
%
%  A + C -> D
%
%where C is the desired product.  The rate of
%the first reaction is given by k1*A*B where
%A and B are the concentrations of the reactants,
%and the rate of the second by k2*A*C.
pause
%We thus get a system of three first order
%differential equations:
%
%  dA/dt = -k1*A*B -k2*A*C
%
%  dB/dt = -k1*A*B
%
%  dC/dt =  k1*A*B -k2*A*C
%
%We want to solve these simultaneously to see
%when the product C will reach maximum concentration.
pause
%We store the differential equations for A, B,
%and C in the function react.m given below.  We
%shall use the matlab ode solver ode23.m to find
%the solution as a function of time.  We pass the
%reaction rate constants as global variables.
global k1 k2
k1=.5;
k2=1;
%We store the concentrations in the array c, such
%that c(1)=A, c(2)=B, and c(3)=C.  We have the
%initial values:
c0=[2,1,0]';
pause
%Thus:
[t,c]=ode23('react',[0,10],c0);
plot(t,c)
xlabel('time')
ylabel('concentration')
legend('A','B','C')
hold on
pause
%We can find the maximum concentration and time by
%using the max command:
[cmax,i]=max(c(:,3))
tmax=t(i)
pause
%We can do a bit better by using the values around
%this maximum to do an interpolation.  We fit a 
%parabola to the maximum value and the values on
%either side of it:
a=[ones(3,1),t(i-1:i+1),t(i-1:i+1).^2];
x=a\c(i-1:i+1,3);
tmax=-x(2)/2/x(3)
cmax=[1,tmax,tmax^2]*x
plot(tmax,cmax,'ow')
%which is a bit higher.  The difference is because
%the ode23 routine evaluates the function only at
%discrete times, thus interpolation is necessary to
%determine the actual maximum.
hold off
echo off




_____________________________



function cprime=react(t,c)
global k1 k2
%This function returns the derivatives of
%the concentrations with respect to time.
%In order to be used with ode23.m it must
%take in time as well as the concentrations,
%although the derivatives are not time
%dependent.
cprime=zeros(size(c));
cprime(1)=-k1*c(1)*c(2) - k2*c(1)*c(3);
cprime(2)=-k1*c(1)*c(2);
cprime(3)=k1*c(1)*c(2) - k2*c(1)*c(3);