% Modified version of the MATLAB template file SolvNonDif.m used in the
% solution of Example 1.4 of A First Course on Kinetics and Reaction
% Engineering.
%
function z = Example_1_4
    % Known quantities and constants (in consistent units)
    n0_CO = 34.;
    n0_H2 = 66.;
    fCO = 0.4;
    R = 1.5;
    S = 8.3;
    
    % Function that evaluates the equations
	function f = evalEqns(z)
        x1 = z(1);
        x2 = z(2);
        x4 = z(3);
		f = [
			fCO*n0_CO - x1 - x4;
            R*(n0_CO - x1 - x4) - (n0_H2 - 2*x1 - 3*x2 - 3*x4);
            x1 + x2 + S*x2;
		];
    end % of internal function evalEqns
    
    % guesses for the solution
    z_guess = [
        0.1
        0.1
        0.1
    ];

    % Solve the set of algebraic equations
    z = fsolve(@evalEqns, z_guess);
    
    f = evalEqns(z)
    
    % Calculate the final composition
    x1 = z(1);
    x2 = z(2);
    x4 = z(3);
    
    nCO = n0_CO - x1 - x4
    nH2 = n0_H2 - 2*x1 - 3*x2 - 3*x4
    nCH3OH = x1 + x2
    nCO2 = -x2
    nH2O = x2 + x4
    nCH4 = x4
    
	
end % of Example_1_4.m