how to I solve 2 eqns with 2 unknowns and plot 5 plots according to these solutions?

clc;

clear;

deltaGrxn_std_rxn1 = -24800; %(J/mol)

deltaHrxn_std_rxn1 = -90200; %(J/mol)

deltaSrxn_std_rxn1 = -219;

deltaGrxn_std_rxn2 = -28600;

deltaHrxn_std_rxn2 = -41200;

deltaSrxn_std_rxn2 = -42;

R = 8.314;

T = 298:10:800;

for j = 1:length(T)

deltaGrxn1_funcT = (deltaHrxn_std_rxn1) - (T(j) .* deltaSrxn_std_rxn1);

Ka_rxn1(j) = exp((-deltaGrxn1_funcT)./(R.*T(j)));

deltaGrxn2_funcT = (deltaHrxn_std_rxn2) - (T(j) .* deltaSrxn_std_rxn2);

Ka_rxn2(j) = exp((-deltaGrxn2_funcT)./(R.*T(j)));

syms ksi1

syms ksi2

y_CO = (1-ksi1-ksi2)./(3-2.*ksi1);

y_H2 = (2-2.*ksi1+ksi2)./(3-2.*ksi1);

y_CH3OH = (ksi1)./(3-2.*ksi1);

y_H2O = (1-ksi2)./(3-2.*ksi1);

y_CO2 = (ksi2)./(3-2.*ksi1);

eqn1 = ((y_CH3OH)./((y_H2O.^2).*(y_CO))) - Ka_rxn1(j) == 0

eqn2 = ((y_CO2.*y_H2)./(y_CO.*y_H2O)) - Ka_rxn2(j) ==0

E = [eqn1, eqn2];

S = vpasolve(E, ksi1, ksi2, [0 1]);

end

plot(T,y_CO)

hold on

plot(T,y_H2)

hold on

plot(T,y_CH3OH)

hold on

plot(T,y_H2O)

hold on

plot(T,y_CO2)

Answer 1

[Alan Stevens]

Here's a possible way. Only you will know if the results make sense

deltaGrxn_std_rxn1 = -24800; %(J/mol)
deltaHrxn_std_rxn1 = -90200; %(J/mol)
deltaSrxn_std_rxn1 = -219; 
deltaGrxn_std_rxn2 = -28600; 
deltaHrxn_std_rxn2 = -41200; 
deltaSrxn_std_rxn2 = -42; 
R = 8.314; 
T = 298:10:800;

deltaGrxn1_funcT = (deltaHrxn_std_rxn1) - (T.* deltaSrxn_std_rxn1);
Ka_rxn1 = exp((-deltaGrxn1_funcT)./(R.*T));
deltaGrxn2_funcT = (deltaHrxn_std_rxn2) - (T.* deltaSrxn_std_rxn2);
Ka_rxn2 = exp((-deltaGrxn2_funcT)./(R.*T));

y_CO = @(ksi1,ksi2) (1-ksi1-ksi2)./(3-2.*ksi1);
y_H2 = @(ksi1,ksi2) (2-2.*ksi1+ksi2)./(3-2.*ksi1);
y_CH3OH = @(ksi1,ksi2) (ksi1)./(3-2.*ksi1);
y_H2O = @(ksi1,ksi2) (1-ksi2)./(3-2.*ksi1);
y_CO2 = @(ksi1,ksi2) ksi2./(3-2.*ksi1);


 yco = zeros(size(T));
 yh2 = zeros(size(T));
 ych3oh = zeros(size(T));
 yh2o = zeros(size(T));
 yco2 = zeros(size(T));
 
 K = [0 1];   % Initial guesses
for j = 1:length(T)
    
    
      K0 = K;  % Update initial guesses   
     
      F = @(K) norm(((y_CH3OH(K(1),K(2)))./((y_H2O(K(1),K(2)).^2).*(y_CO(K(1),K(2))))) - Ka_rxn1(j))...
         + norm(((y_CO2(K(1),K(2)).*y_H2(K(1),K(2)))./(y_CO(K(1),K(2)).*y_H2O(K(1),K(2)))) - Ka_rxn2(j));

      K = fminsearch(F,K0);
      yco(j) = y_CO(K(1),K(2));
      yh2(j) = y_H2(K(1),K(2));
      ych3oh(j) = y_CH3OH(K(1),K(2));
      yh2o(j) = y_H2O(K(1),K(2));
      yco2(j) = y_CO2(K(1),K(2));
end
    
    figure(1)
    plot(T,yh2,T,yco2),grid
    xlabel('T'), ylabel('y')
    legend('H2','CO2')

    
    figure(2)
    plot(T,yco,T,yh2o),grid
    xlabel('T'),ylabel('y')
    legend('CO','H2O')

    
    figure(3)
    plot(T,ych3oh),grid
    xlabel('T'),ylabel('y_CH3OH')

CO and H2O are numerically too close to separate graphically.