function Example_38_1
%Example_38_1 Code used in the solution of Example 38.1 of AFCoKaRE

    % Known quantities and constants
    k = 2.0E-2; % /s
    eps = 0.4;
    VFR = 1.0; % cm^3/s
    CA0 = 1.16E-3; % mol/cm^3
    Rp = 0.3/2; % cm
    DA = 2.7E-7; % cm^2/s
    fA = 0.85;
    D = 2.5; % cm
    
    % Calculate the effectiveness factor
    Thiele = Rp*sqrt(k/DA) % solution, equation 5
    eff_factor = 3/Thiele*(1/tanh(Thiele) - 1/Thiele) % solution, equation 6
    %eff_factor = 1
    
	% Function that evaluates the ODEs
    function dzdt = odeqns(t,z)
        dzdt = [
            -pi()*D^2/4*k*eff_factor*(1-eps)*z(1)/VFR; 
        ]; % solution, equations 1 and 4
    end % of internal function odeqns
        
    % Initial values
    t0 = 0;
    z0 = [
        CA0*VFR % solution, equation 2
    ];

    tf = 460.; 
    %tf = 32.2
    options = odeset('Events',@stop);
	[t, zz, te, ze, ie] = ode45(@odeqns,[t0, tf],z0,options);
    
    % Function that provides the integration stopping criterion
    function [stop_when, isterminal, direction] = stop(t,z)
        isterminal = 1;
        direction = 0;
        
        stop_when = (1-fA)*CA0*VFR - z(1); % solution, equation 7
    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. Set tf to a larger value and re-execute');
        display('    the function')
    end % if (t_f >= tf)
	
    check = (CA0*VFR - z(1))/CA0/VFR
    L = t_f

end % of Example_38_1
