echo off
clc
% 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 Adams
% Method though one sees the BDF Method takes significantly fewer
% steps.  Even though a BDF step is more expensive due to the Newton
% iteration to "solve" the implicit equation, the overall performance
% of the BDF surpasses that of the Adams.
%
% 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.
%
c Chemistry is described by:
c
c   Cl + H2  -->  HCl + H       k1 = 1.6e-14
c   H + Cl2  -->  HCl + Cl      k2 = 2.0e-11
c   H + O2   -->  HO2           k3 = 3.6e-13
c   Cl + O2  -->  ClO2          k4 = 1.3e-14
c   Cl + ClO2-->  Cl2 + O2      k5 = 1.4e-10
c
c Labelling the concentrations for the dependent variable:
c
c   Y(1) = Cl
c   Y(2) = H2
c   Y(3) = Cl2
c   Y(4) = HCl
c   Y(5) = H
c   Y(6) = O2
c   Y(7) = HO2
c   Y(8) = ClO2
c
c We have the rate differential equations (from the subroutine ozone):
c
c      ydot(1) = -k(1)*y(1)*y(2)
c     1          +k(2)*y(3)*y(5)
c     2          -k(4)*y(1)*y(6)
c     3          -k(5)*y(1)*y(8)
c      ydot(2) = -k(1)*y(1)*y(2)
c      ydot(3) = -k(2)*y(3)*y(5)
c     1          +k(5)*y(1)*y(8)
c      ydot(4) =  k(1)*y(1)*y(2)
c     1          +k(2)*y(3)*y(5)
c      ydot(5) =  k(1)*y(1)*y(2)
c     1          -k(2)*y(3)*y(5)
c     2          -k(3)*y(5)*y(6)
c      ydot(6) = -k(3)*y(5)*y(6)
c     1          -k(4)*y(1)*y(6)
c     2          +k(5)*y(1)*y(8)
c      ydot(7) =  k(3)*y(5)*y(6)
c      ydot(8) =  k(4)*y(1)*y(6)
c     1          -k(5)*y(1)*y(8)

% pass the reaction coefficients to the derivative module
$ ozone.m
k1=1.6e-14; k2=2.0e-11; k3=3.6e-13; k4=1.3e-14; k5=1.4e-10;
global k1 k2 k3 k4 k5;
%

% don't need since canode23 has a break tfinal = input('enter final time');
tfinal = 100;

z0 = [0 0 50 angle];
t0 = 0;
tol = 1.e-7;
trace = 0;
[t,z] = canode23('cannon',t0,tfinal,z0,tol,trace);
%
mntx=[50 90 130]; mnty=[0 60 0];
plot(z(:,1),z(:,2),mntx,mnty),
stit=['Cannon shot    Launch angle=' num2str(angled) ' degrees;  Accuracy='num2str(tol)];
title(stit),xlabel('Distance'),ylabel('Height')