% MATLAB template file for solving a set of coupled initial value ordinary
% differential equations of the form
%   dz1/dt = f1(t, z1, z2, z3, ..., zn); z1(0) = c1
%   dz2/dt = f2(t, z1, z2, z3, ..., zn); z2(0) = c2
%              .
%              .
%              .
%   dzn/dt = fn(t, z1, z2, z3, ..., zn); zn(0) = cn
% where c1 through cn are known constants and the final value of one
% dependent variable, zi, is known
%
% See Supplemental Unit S8 of "A First Course on Kinetics and Reaction
% Engineering" and/or the MATLAB documentation for additional information
% or if questions arise.
%
% Carl R. F. Lund, 4/22/13
%   last rev. 6/19/13
%
% Returned quantities
%   t_f (scalar) - the value of the independent variable at which the
%       values of the dependent variables were determined
%   z (column vector) - values of z1, z2, ..., zn found by integrating the
%       differential equations from the specified initial conditions to the
%       value of the independent variable returned as t_f
%
% This file requires six modifications before each use; the locations
% where editing is required are indicated by the comment "% EDIT HERE"
%
function [t_f,z] = SolvIVDifD
    % Known quantities and constants (in consistent units)
% EDIT HERE (Required modification 1 of 6): enter known quantities
    

	% Function that evaluates the ODEs
    function dzdt = odeqns(t,z)
% EDIT HERE (Required modification 2 of 6):
        % t is the value of the independent variable
        % z is a vector containing the values of the dependent variables
        dzdt = [
            % Evaluate dz1/dt = f1(t, z1, z2, z3, ..., zn) here
			% Evaluate dz2/dt = f2(t, z1, z2, z3, ..., zn) here
            % and so on through fn, one per line
        ];
    end % of internal function odeqns
        
    % Initial values
% EDIT HERE (Required modification 3 of 6):
    t0 = ;% insert the initial value of independent variable here
    z0 = [
        % insert the initial values of dependent variables z1 here
        % insert initial values for z2, z3, ..., one per line
    ];

% EDIT HERE (Required modification 4 of 6):
    % Set the final value, tf, of the independent variable on the next
    % line to a value that is sufficiently large to cause the known final
    % dependent variable value to be reached before the independent
    % variable reaches tf.
    
    tf = ;% insert a large number, as just instructed
    options = odeset('Events',@stop);
	[t, zz, te, ze, ie] = ode45(@odeqns,[t0, tf],z0,options);
    % NOTE: For stiff equations replace ode45 with ode15s
    
    % Function that provides the integration stopping criterion
    function [stop_when, isterminal, direction] = stop(t,z)
        isterminal = 1;
        direction = 0;
        
% EDIT HERE (Required modification 5 of 6) The variable stop_when should 
% equal zero when the desired stopping criterion is reached
        stop_when = ;% replace this comment with the stopping criterion
    end % of internal function stop

    % Extract the final values of the independent and dependent variables
	last_value = length(t);
	t_f = t(last_value); % final value of the independent variable
	z = zz(last_value,:); % final values of the dependent variables
    
    % Warn if the independent variable reached tf above
    if (t_f >= tf)
        display('WARNING: the final value of the independent variable was')
        display('reached, instead of the final value of the dependent')
        display('variable.');
        display('Set tf to a larger value and re-execute the function')
    end
	
% EDIT HERE (Required modification 6 of 6):
    % Enter code to calculate any other desired quantities using the 
    % results contained in z. (Alternatively, z will be returned when this
    % template file terminates; other quantities can then be calculated
    % at the MATLAB command prompt using the returned values
    
end % of SolvIVDifD.m
