% Computes the equilibrium composition of a reacting mixture

  global ng ns nr tempin prin ctot xg nu
     global keq298 keq delhr 
     global n ntot 
    global rgas
%number of species and number of reactions
ng = 4 ;
nr = 2 ;
% initialize the vectors.
keq298 = zeros(nr, 1 );
keq = zeros(nr, 1 );
delhr = zeros(nr, 1 );
xg = zeros(ng,1);
nu = zeros ( nr, ng);
zeta = zeros(nr,1) ;
% provide values for stoichiometry and eq constants.
nu(1,1) = -1.0 ;
nu(1,2) = -1.0;
nu(1,3) = 2.0 ;
nu(1,4)= 0.0;
nu(2,1) = 0.0 ;
nu(2,2) = -0.5;
nu(2,3) = -1.0 ;
nu(2,4)= 1.0;
keq298(1) = 4.069E-31;
delhr(1) = 180500.0; 
keq298(2) = 1.22E+06 ;
delhr(2) = -57070.0;
 %assumed constant (not a function of T)
rgas = 8.314472 ; % gas constant; use consistent units
% provide the feed conditions
tempin = 1800.0;
prin = 1.0e+05;
% initial moles
xg(1) = 96.7;
xg(2) = 3.3;
xg(3) = 0.0;
xg(4) =0.0; 
% provide initial values for extent of reactions
zeta(1) = 0.4;
zeta(2) = 0.1;
zeta = fsolve ('fun_eq', zeta)
----------------------------------------------------
function fvec = fun_eq( zeta)
% Defines the (nr) functions to be solved
          global ng ns nr tempin prin ctot xg nu
     global keq298 keq delhr 
     global n ntot 
    global rgas
fvec = zeros (nr,1) ;
ptotatm = prin /1.0e+05 ;

% find K at the desired conditions
for i = 1: nr
keq(i)=keq298(i)*exp( delhr(i) /rgas *(1./298-1./tempin)) ;
end
%number of moles of each species for the given extents
% nj0 = 1 for these calculations
for j = 1: ng
n(j) = xg(j);
for i=1:nr
n(j) = n(j) + nu(i,j) * zeta(i);
end
end
% find ntot
ntot = 0.0;
for j = 1: ng
ntot = ntot + n(j);
end
% find the mole fractions for each species
pp = n /ntot * ptotatm;
% set up discrepancy function for each reaction
for i = 1: nr
prod = 1.0;
for j = 1: ng
prod = prod * pp(j)^nu(i,j);
end
fvec(i) = keq(i) - prod
end