% Modified version of the MATLAB template file SolvNonDif.m used to solve
% Example 3.2 of "A First Course on Kinetics and Reaction Engineering."
%
function Example_3_2
    % Known quantities and constants (in mol atm)
    % stoichiometric coefficients
    nu_co = -1;
    nu_h2 = -2;
    nu_meoh = 1;
    % initial mole fractions
    y0_co = 0.25;
    y0_h2 = 0.5;
    y0_meoh = 0.0;
    y0_co2 = 0.25;
    P = 5.0/.10133; % converted to atm
    K_500 = 0.006153 ;       
    n0_tot = 1.0; % basis
    
    % Calculate the initial moles
    n0_co =  y0_co*n0_tot;
    n0_h2 =  y0_h2*n0_tot;
    n0_meoh =  y0_meoh*n0_tot;
    n0_co2 =  y0_co2*n0_tot;
    
    
    % Function that evaluates the equations
	function f = evalEqns(z)
        extent = z(1);
        n_co = n0_co - extent;
        n_h2 = n0_h2 - 2*extent;
        n_meoh = n0_meoh + extent;
        n_co2 = n0_co2;
        n_tot = n_co + n_h2 + n_meoh + n_co2;
        y_co = n_co/n_tot;
        y_h2 = n_h2/n_tot;
        y_meoh = n_meoh/n_tot;
		f = [
			K_500*P^2*y_co*y_h2^2 - y_meoh;
		];
    end % of internal function evalEqns
    
    % guesses for the solution

    z_guess = [
        0.01
    ];

    % Solve the set of algebraic equations
    z = fsolve(@evalEqns, z_guess);
    
    % Calculate the final mole fractions
    extent = z(1);
    n_co = n0_co - extent;
    n_h2 = n0_h2 - 2*extent;
    n_meoh = n0_meoh + extent;
    n_co2 = n0_co2;
    n_tot = n_co + n_h2 + n_meoh + n_co2;
    y_co = n_co/n_tot
    y_h2 = n_h2/n_tot
    y_meoh = n_meoh/n_tot
    y_co2 = n_co2/n_tot
	
end % of Example_3_2