function Example6B
clear, clc, format short g, format compact
xguess = [0.1 0.1 0.1]; % initial guess vector 
disp('Variable values at the initial estimate');
fguess=MNLEfun(xguess,330);
disp(' Variable        Value     Function Value')
for i=1:size(xguess,2);
disp([' x' int2str(i) '                 ' num2str(xguess(i)) '            ' num2str(fguess(i))]);
end
options = optimset('Diagnostics',['off'],'TolFun',[1e-9],'TolX',[1e-9]);
CA0 = 1.5; 
CB0 = 1.5; 
for i=1:41 
    T(i,1)=330+(i-1);
    xsolv=fsolve(@MNLEfun,xguess,options,T(i));
    xguess=xsolv;
    CD(i,1)=xsolv(1); CX(i,1)=xsolv(2); CZ(i,1)=xsolv(3);
    CY (i,1)= CX (i)+ CZ(i); 
    CC(i,1) = CD(i) - CY(i); 
    CA(i,1) = CA0 - CD(i) - CZ(i); 
    CB(i,1) = CB0 - CD(i) - CY(i); 
end
disp('                 Complex Chemical Equilibrium, Reactants Concentrations');
disp('                                   Tabular Results');
disp('');
disp('        T(K)          CA         CB  ');
disp([T CA CB ]);
pause
disp('');
disp('                 Complex Chemical Equilibrium, Products Concentrations');
disp('                                   Tabular Results');
disp('');
disp('        T(K)           CC           CD            CX              CY            CZ');
disp([T CC CD CX CY CZ]);;
pause
plot(T,CA,'-',T,CB,'+');
legend(' CA',' CB');
title(' Reactants Equilibrium Concentration')
xlabel('Temperature (K)');
ylabel('Concentration (g-mol/L)');
pause
plot(T,CC,'-',T,CD,'+', T,CX,'*',T,CY,'x', T,CZ,'o');
legend(' CC',' CD',' CX',' CY',' CZ');
title(' Products Equilibrium Concentration')
xlabel('Temperature (K)');
ylabel('Concentration (g-mol/L)');
pause

%- - - - - - - - - - - - - - - - - - - - - -
function fx = MNLEfun(x,T); 
CD = x(1); 
CX = x(2); 
CZ = x(3); 
%T = 330; 
KC1 = 1 / exp(-2396.301 * (1 / 330 - (1 / T))) * .7; 
KC2 = 1 / exp(2421.518 * (1 / 330 - (1 / T))) * 4; 
KC3 = 1 / exp(-8954.7 * (1 / 330 - (1 / T))) * 1 * T / 330; 
CY = CX + CZ; 
CA0 = 1.5; 
CB0 = 1.5; 
CC = CD - CY; 
CA = CA0 - CD - CZ; 
CB = CB0 - CD - CY; 
fx(1,1) = CC * CD - (KC1 * CA * CB); 
fx(2,1) = CX * CY - (KC2 * CB * CC); 
fx(3,1) = CZ - (KC3 * CA * CX);