echo on
% This problem is stated in p. 337 of "Numerical Methods and Software"
% by Kahaner, Moler and Nash; Prentice-Hall, 1989 (p8-14)
%
% It is "mildly stiff", therefore it is still solvable by the Runge-Kutta
%
% The reaction time constants for this problem (below) only vary by 4
% orders of magnitude.  Stiffer problems typically have a much wider
% variation in reaction time constants.
%
% Chemistry is described by:
%
%   Cl + H2  -->  HCl + H       k1 = 1.6e-14
%   H + Cl2  -->  HCl + Cl      k2 = 2.0e-11
%   H + O2   -->  HO2           k3 = 3.6e-13
%   Cl + O2  -->  ClO2          k4 = 1.3e-14
%   Cl + ClO2-->  Cl2 + O2      k5 = 1.4e-10
%
% Labelling the concentrations for the dependent variable:
%
%   Y(1) = Cl
%   Y(2) = H2
%   Y(3) = Cl2
%   Y(4) = HCl
%   Y(5) = H
%   Y(6) = O2
%   Y(7) = HO2
%   Y(8) = ClO2
%
% We have the rate differential equations (from the subroutine ozone):
%
%      ydot(1) = -k(1)*y(1)*y(2)
%     1          +k(2)*y(3)*y(5)
%     2          -k(4)*y(1)*y(6)
%     3          -k(5)*y(1)*y(8)
%      ydot(2) = -k(1)*y(1)*y(2)
%      ydot(3) = -k(2)*y(3)*y(5)
%     1          +k(5)*y(1)*y(8)
%      ydot(4) =  k(1)*y(1)*y(2)
%     1          +k(2)*y(3)*y(5)
%      ydot(5) =  k(1)*y(1)*y(2)
%     1          -k(2)*y(3)*y(5)
%     2          -k(3)*y(5)*y(6)
%      ydot(6) = -k(3)*y(5)*y(6)
%     1          -k(4)*y(1)*y(6)
%     2          +k(5)*y(1)*y(8)
%      ydot(7) =  k(3)*y(5)*y(6)
%      ydot(8) =  k(4)*y(1)*y(6)
%     1          -k(5)*y(1)*y(8)

% type ozone
% pause % Strike any key to continue.
% clc
% clear
%	To simulate the differential equation defined in ozone over the
%	interval  0 < t < ?, we invoke ODE23:
%
k1=1.6e-14; k2=2.e-11; k3=3.6e-13; k4=1.3e-14; 1.4e-10;
global k1 k2 k3 k4 k5;

t0 = 0;  tfinal=0.0001;
y0 = [1.E14 1.62E18 3.25E16 0 0 4.84E18 0 0]';
	
% Define initial conditions.

tol = 1.e-4;            % Accuracy
trace = 0; % turn off trace

[t,y] = ode23('ozone',t0,tfinal,y0,tol,trace);

plot(t,y(:,1),t,y(:,7),t,y(:,8)), ...
title('ozone chemistry equation time history'), pause
clc