% Modified version of the AFCoKaRE template file SolvIVDifD.m used to solve
% Activity 26.1.
%
function Activity_26_1
    % Known quantities and constants (in consistent units)
    yA0 = 0.05;
    yB0 = 0.2;
    yI0 = 0.75;
    T0 = 350 + 273.15;
    P0 = 2;
    fA = 0.8;
    k0 = 1.4e6;
    E = 30;
    dH = -25;
    Cp = 1.987e-3;
    nY0 = 0;
    nZ0 = 0;
    ntot0 = 100; %Basis
    nA0 = yA0*ntot0;
    nB0 = yB0*ntot0;
    nI0 = yI0*ntot0;
    nI = nI0;
    P = P0;
    Rpv = 0.08206;
    Ren = 1.987e-3;
    VFR0 = ntot0*Rpv*T0/P0;
    nAf = nA0*(1-fA);

	% Function that evaluates the ODEs
    function dzdt = odeqns(t,z)
        nA = z(1);
        nB = z(2);
        nY = z(3);
        nZ = z(4);
        T = z(5);
        nTot = nA + nB + nY + nZ + nI;
        CA = nA*P/nTot*Rpv*T;
        r1 = k0*exp(-E/Ren/T)*CA;
        dzdt = [
            -2*r1;
            -r1;
            2*r1;
            2*r1;
            -r1*dH/nTot/Cp;
        ];
    end % of internal function odeqns
        
    % Initial values
    t0 = 0;
    z0 = [
        nA0;
        nB0;
        nY0;
        nZ0;
        T0;
    ];

    tf = 1e10;
    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;
        % final condition
        stop_when = z(1)-nAf;
    end % of internal function stop

    % Return the values of the independent and dependent variables when the
    % integration stopping condition was satisfied
	last_value = length(t);
	t_f = t(last_value);
	z = zz(last_value,:);
	
    T_f_K = z(5)
    tau_min = t_f/VFR0
end % of Activity_26_1.m
