% Modified version of the MATLAB template file SolvIVDifI.m used in the
% solution of Activity 19.1 of "A First Course on Kinetics and Reaction
% Engineering."
%
function Activity_19_1(t_final)
    % Known quantities and constants (L, K, atm, s, kJ, mol)
    ppmA0 = 1500;
    ppmB0 = 1000;
    pctC0 = 7;
    V = 3;
    T0 = 1115;
    P0 = 1.7;
    k01 = 6.1e16;
    k02 = 5.5e13;
    E1 = 250;
    E2 = 320;
    dH1 = -1700;
    dH2 = -800;
    cpI = 32/1000;
    R_pv = 0.08206;
    R_kJ = 8.31446/1000;
    % Initial moles
    n0 = P0*V/R_pv/T0;
    nA0 = ppmA0*n0/1e6;
    nB0 = ppmB0*n0/1e6;
    nC0 = pctC0*n0/100;
    nI = n0-nA0-nB0-nC0;
    
	% Function that evaluates the ODEs
    function dzdt = odeqns(t,z)
        nA = z(1);
        nB = z(2);
        nC = z(3);
        nY = z(4);
        nZ = z(5);
        T = z(6);
        r1 = k01*exp(-E1/R_kJ/T)*nA/V*nB/V;
        r2 = k02*exp(-E2/R_kJ/T)*nA/V;
        dzdt = [
            V*(-4*r1 - 4*r2)
            V*(-4*r1 + 4*r2)
            V*(-r1 - 5*r2)
            V*4*r1
            V*(6*r1 + 6*r2)
            -V*(r1*dH1 + r2*dH2)
            0
        ];
    end % of internal function odeqns

    % function that evaluates the mass matrix
    function M = mass(t,z)
        M = zeros(7,7);
        M(1,1) = 1;
        M(2,2) = 1;
        M(3,3) = 1;
        M(4,4) = 1;
        M(5,5) = 1;
        M(6,6) = (z(1)+z(2)+z(3)+z(4)+z(5)+nI)*cpI;
        M(6,7) = -V;
        M(7,1) = R_kJ*z(6);
        M(7,2) = R_kJ*z(6);
        M(7,3) = R_kJ*z(6);
        M(7,4) = R_kJ*z(6);
        M(7,5) = R_kJ*z(6);
        M(7,6) = R_kJ*(z(1)+z(2)+z(3)+z(4)+z(5)+nI);
        M(7,7) = -V;
    end % of internal function mass
    
        
    % Initial and final values
    t0 = 0;
    z0 = [
        nA0
        nB0
        nC0
        0
        0
        T0
        P0
    ];
    tf = t_final;
    
    % Solve the ODEs
    options = odeset('Mass',@mass);
    [t,zz] = ode45(@odeqns,[t0 tf],z0,options);
% NOTE: For stiff equations replace ode45 with ode15s

    % 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,:);
	
    % Calculate ppm of B
    ppmB = 1e6*z(2)/(z(1)+z(2)+z(3)+z(4)+z(5)+nI)
    
end % of Activity_19_1.m
