RCM Computation using UNIQUAC Model
Case of a Quaternary Reactive System : 
Acetic Acid - Butanol - Water- Butyl Acetate
Author's Data: Housam BINOUS
Department of Chemical Engineering
National Institute of Applied Sciences and Technology
Tunis, TUNISIA
Email: binoushousam@yahoo.com

Off[General::"spell1"]

Vapor pressure and normal boiling points

Psat1 = 
      10^3  Exp[A1 + B1/(C1 + T) + D1 Log[T] + E1 T^F1]

Psat2 = 
      10^3 Exp[A2 + B2/(C2 + T) + D2 Log[T] + E2 T^F2]

Psat3 = 
      10^3 Exp[A3 + B3/(C3 + T) + D3 Log[T] + E3 T^F3]

Psat4 = 
      10^3  Exp[A4 + B4/(C4 + T) + D4 Log[T] + E4 T^F4]

A1 = 61.34090; B1 = -6768.88; C1 = 0; D1 = -6.72663; 
    E1 = 4.84 10^-6; F1 = 2;

A2 = 108.8; B2 = -1.007 10^4; C2 = 0; D2 = -13.26; 
    E2 = 4.38 10^-6; F2 = 2;

A3 = 103.8; B3 = -8764; C3 = 0; D3 = -13.12; 
    E3 = 8.44 10^-6; F3 = 2;

A4 = 65.92780; B4 = -7227.53; C4 = 0; D4 = -7.17695; 
    E4 = 4.03 10^-6; F4 = 2;

FindRoot[Psat1 ==  1.01325 10^5, {T, 380}]

390.951 - 273.15

FindRoot[Psat2 ==  1.01325 10^5, {T, 350}]

392.37 - 273.15

FindRoot[Psat3 ==  1.01325 10^5, {T, 380}]

399.512 - 273.15

FindRoot[Psat4 ==  1.01325 10^5, {T, 380}]

373.153 - 273.15

Activity coefficients using UNIQUAC model
Component 1 = Acetic Acid
Component 2 = Butanol
Component 3 = Butyl Acetate
Component 4 = Water

tau[j_, i_] := Exp[-(U[j, i] - U[i, i])/(R T)]

U[1, 1] = 0.0; U[1, 2] = -157.20; U[1, 3] = -298.434204; 
    U[1, 4] = 427.740509;

U[2, 1] = 177.1036; U[2, 2] = 0.0; U[2, 3] = 279.02160; 
    U[2, 4] = -36.265091;

U[3, 1] = 712.23492 ; U[3, 2] = -111.1324; U[3, 3] = 0.0; 
    U[3, 4] = 738.16931;

U[4, 1] = -305.451507; U[4, 2] = 727.3850; U[4, 3] = 398.29040; 
    U[4, 4] = 0.0;

Gam1 = Exp[Log[Gam1C] + Log[Gam1R]]

Gam2 = Exp[Log[Gam2C] + Log[Gam2R]]

Gam3 = Exp[Log[Gam3C] + Log[Gam3R]]

Gam4 = Exp[Log[Gam4C] + Log[Gam4R]]

Gam1C = 
        Exp[Log[phi[1]/x[1]] + z/2 q[1] Log[v[1]/phi[1]] + l[1] - 
            phi[1]/x[1] (∑_(j = 1)^4) x[j] l[j]];

Gam2C = 
        Exp[Log[phi[2]/x[2]] + z/2 q[2] Log[v[2]/phi[2]] + l[2] - 
            phi[2]/x[2] (∑_(j = 1)^4) x[j] l[j]];

Gam3C = 
        Exp[Log[phi[3]/x[3]] + z/2 q[3] Log[v[3]/phi[3]] + l[3] - 
            phi[3]/x[3] (∑_(j = 1)^4) x[j] l[j]];

Gam4C = 
        Exp[Log[phi[4]/x[4]] + z/2 q[4] Log[v[4]/phi[4]] + l[4] - 
            phi[4]/x[4] (∑_(j = 1)^4) x[j] l[j]];

r[1] = 2.20230; r[2] = 3.45420; r[3] = 4.82740; 
    r[4] = 0.920000;

q[1] = 2.07200; q[2] = 3.05200; q[3] = 4.19600; 
    q[4] = 1.39970;

l[i_] := z/2 (r[i] - q[i]) - (r[i] - 1)

z = 10;

Gam1R = 
        Exp[q[1] (1 - 
                Log[(∑_(j = 1)^4) v[j] tau[j, 
                        1]] - (∑_(j = 1)^4)(( v[j] tau[1, j])/((∑_(k = 1)^4) v[k] tau[k, j])))];

Gam2R = 
        Exp[q[2] (1 - 
                Log[(∑_(j = 1)^4) v[j] tau[j, 
                        2]] - (∑_(j = 1)^4)(( v[j] tau[2, j])/((∑_(k = 1)^4) v[k] tau[k, j])))];

Gam3R = 
        Exp[q[3] (1 - 
                Log[(∑_(j = 1)^4) v[j] tau[j, 
                        3]] - (∑_(j = 1)^4)(( v[j] tau[3, j])/((∑_(k = 1)^4) v[k] tau[k, j])))];

Gam4R = 
        Exp[q[4] (1 - 
                Log[(∑_(j = 1)^4) v[j] tau[j, 
                        4]] - (∑_(j = 1)^4)(( v[j] tau[4, j])/((∑_(k = 1)^4) v[k] tau[k, j])))];

phi[i_] := r[i] x[i]/(∑_(j = 1)^4) r[j] x[j]

v[i_] := q[i] x[i]/(∑_(j = 1)^4) q[j] x[j]

R = 1.98721;

Dimerization constant for acetic acid

K = 10^(-12.5454 + 3166/(T));

Total pressure

P = 1.01325 10^5;

Chemical reaction equilibrium constant

Keq = 12.5;

Computing the residue curve map. This system exhibits a reactive azeotrope.

t = 0; h = 0.5; X[1, 0] = 0.7; X[2, 0] = 0.7; x[1, -1] = 0.1; 
    x[2, -1] = 0.2; x[3, -1] = 0.1; x[4, -1] = 0.8; 
    y[1, -1] = 0.5;

While[
        t < 500, 
{
sol = 
            Chop[FindRoot[{X[1, t] ==  
                    x[1] + x[3], 
X[2, t] ==  
                    x[2] + x[3], 
                  Gam1 Gam2 x[1] x[2] - 
                      Gam3 Gam4 x[3] x[4]/Keq ==  
                    0, 
x[1] + x[2] + x[3] + 
                      x[4] ==  
                    1, 
Psat1 Gam1 x[
                            1]/((1 + 
                                  Sqrt[(1 + 4 K Psat1)])/(1 + 
                                  Sqrt[(1 + 
                                        4 K P y[
                                        1] (2 - y[1]))])) + 
                      Psat2 Gam2 x[
                            2]/((2 (1 - y[1] + 
                                      Sqrt[(1 + 
                                        4 K P y[
                                        1] (2 - 
                                        y[1]))]))/((2 - 
                                      y[1]) (1 + 
                                      Sqrt[(1 + 
                                        4 K P y[
                                        1] (2 - y[1]))]))) + 
                      Psat3 Gam3 x[
                            3]/((2 (1 - y[1] + 
                                      Sqrt[(1 + 
                                        4 K P y[
                                        1] (2 - 
                                        y[1]))]))/((2 - 
                                      y[1]) (1 + 
                                      Sqrt[(1 + 
                                        4 K P y[
                                        1] (2 - y[1]))]))) + 
                      Psat4 Gam4 x[
                            4]/((2 (1 - y[1] + 
                                      Sqrt[(1 + 
                                        4 K P y[
                                        1] (2 - 
                                        y[1]))]))/((2 - 
                                      y[1]) (1 + 
                                      Sqrt[(1 + 
                                        4 K P y[
                                        1] (2 - 
                                        y[1]))]))) ==  
                    P, 
y[1] ==  
                    Psat1 Gam1 x[
                          1]/((1 + Sqrt[(1 + 4 K Psat1)])/(1 + 
                                Sqrt[(1 + 
                                      4 K P y[
                                        1] (2 - y[1]))])) 1/
                        P}, {T, 350}, {x[3], x[3, t - 1]}, {x[4], 
                  x[4, t - 1]}, {x[1], x[1, t - 1]}, {x[2], 
                  x[2, t - 1]}, {y[1], y[1, t - 1]}, 
                MaxIterations -> 100000]], 
T[t] = 
            Chop[sol[[1, 2]]], 
x[3, t] = 
            Chop[sol[[2, 2]]], 
          x[4, t] = Chop[sol[[3, 2]]], 
x[1, t] = 
            Chop[sol[[4, 2]]], 
x[2, t] = 
            Chop[sol[[5, 2]]], 
y[1, t] = 
            Chop[sol[[6, 2]]], 
y[2, 
              t] = Chop[
                        Psat2 Gam2 x[2, 
                              t]/((2 (1 - y[1, t] + 
                                        Sqrt[(1 + 
                                        4 K P y[1, 
                                        t] (2 - y[1, 
                                        t]))]))/((2 - 
                                        y[1, t]) (1 + 
                                        Sqrt[(1 + 
                                        4 K P y[1, 
                                        t] (2 - y[1, t]))]))) 1/
                            P] /. x[3] -> x[3, t] /. x[4] -> x[4, t] /. 
                  x[1] -> x[1, t] /. x[2] -> x[2, t] /. 
              T -> T[t], 
y[3, 
              t] = Chop[
                        Psat3 Gam3 x[3, 
                              t]/((2 (1 - y[1, t] + 
                                        Sqrt[(1 + 
                                        4 K P y[1, 
                                        t] (2 - y[1, 
                                        t]))]))/((2 - 
                                        y[1, t]) (1 + 
                                        Sqrt[(1 + 
                                        4 K P y[1, 
                                        t] (2 - y[1, t]))]))) 
                          1/P] /. x[3] -> x[3, t] /. x[4] -> x[4, t] /. 
                  x[1] -> x[1, t] /. x[2] -> x[2, t] /. 
              T -> T[t], 
Y[1, t] = 
            Chop[y[1, t] + y[3, t]], 
Y[2, t] = 
            Chop[y[2, t] + y[3, t]], 
X[1, t + 1] = 
            Chop[X[1, t] - h (X[1, t] - Y[1, t])], 
X[
              2, t + 1] = Chop[X[2, t] - h (X[2, t] - Y[2, t])], 
          Print[t], t++}] // Timing

tble1 = Table[{X[1, t], X[2, t]}, {t, 0, 499}];

plt1 = 
      ListPlot[tble1, PlotRange -> {{0, 1}, {0, 1}}, 
        PlotStyle -> RGBColor[1, 0, 0], PlotJoined ->  True, 
        Frame -> True]

t = 0; h = 0.005; X[1, 0] = 0.7; X[2, 0] = 0.7; x[1, -1] = 0.1; 
    x[2, -1] = 0.2; x[3, -1] = 0.1; x[4, -1] = 0.8; 
    y[1, -1] = 0.5;

While[
        t < 500, 
{
sol = 
            Chop[FindRoot[{X[1, t] ==  
                    x[1] + x[3], 
X[2, t] ==  
                    x[2] + x[3], 
                  Gam1 Gam2 x[1] x[2] - 
                      Gam3 Gam4 x[3] x[4]/Keq ==  
                    0, 
x[1] + x[2] + x[3] + 
                      x[4] ==  
                    1, 
Psat1 Gam1 x[
                            1]/((1 + 
                                  Sqrt[(1 + 4 K Psat1)])/(1 + 
                                  Sqrt[(1 + 
                                        4 K P y[
                                        1] (2 - y[1]))])) + 
                      Psat2 Gam2 x[
                            2]/((2 (1 - y[1] + 
                                      Sqrt[(1 + 
                                        4 K P y[
                                        1] (2 - 
                                        y[1]))]))/((2 - 
                                      y[1]) (1 + 
                                      Sqrt[(1 + 
                                        4 K P y[
                                        1] (2 - y[1]))]))) + 
                      Psat3 Gam3 x[
                            3]/((2 (1 - y[1] + 
                                      Sqrt[(1 + 
                                        4 K P y[
                                        1] (2 - 
                                        y[1]))]))/((2 - 
                                      y[1]) (1 + 
                                      Sqrt[(1 + 
                                        4 K P y[
                                        1] (2 - y[1]))]))) + 
                      Psat4 Gam4 x[
                            4]/((2 (1 - y[1] + 
                                      Sqrt[(1 + 
                                        4 K P y[
                                        1] (2 - 
                                        y[1]))]))/((2 - 
                                      y[1]) (1 + 
                                      Sqrt[(1 + 
                                        4 K P y[
                                        1] (2 - 
                                        y[1]))]))) ==  
                    P, 
y[1] ==  
                    Psat1 Gam1 x[
                          1]/((1 + Sqrt[(1 + 4 K Psat1)])/(1 + 
                                Sqrt[(1 + 
                                      4 K P y[
                                        1] (2 - y[1]))])) 1/
                        P}, {T, 350}, {x[3], x[3, t - 1]}, {x[4], 
                  x[4, t - 1]}, {x[1], x[1, t - 1]}, {x[2], 
                  x[2, t - 1]}, {y[1], y[1, t - 1]}, 
                MaxIterations -> 1000000, 
                AccuracyGoal ->  2]], 
T[t] = 
            Chop[sol[[1, 2]]], 
x[3, t] = 
            Chop[sol[[2, 2]]], 
          x[4, t] = Chop[sol[[3, 2]]], 
x[1, t] = 
            Chop[sol[[4, 2]]], 
x[2, t] = 
            Chop[sol[[5, 2]]], 
y[1, t] = 
            Chop[sol[[6, 2]]], 
y[2, 
              t] = Chop[
                        Psat2 Gam2 x[2, 
                              t]/((2 (1 - y[1, t] + 
                                        Sqrt[(1 + 
                                        4 K P y[1, 
                                        t] (2 - y[1, 
                                        t]))]))/((2 - 
                                        y[1, t]) (1 + 
                                        Sqrt[(1 + 
                                        4 K P y[1, 
                                        t] (2 - y[1, t]))]))) 1/
                            P] /. x[3] -> x[3, t] /. x[4] -> x[4, t] /. 
                  x[1] -> x[1, t] /. x[2] -> x[2, t] /. 
              T -> T[t], 
y[3, 
              t] = Chop[
                        Psat3 Gam3 x[3, 
                              t]/((2 (1 - y[1, t] + 
                                        Sqrt[(1 + 
                                        4 K P y[1, 
                                        t] (2 - y[1, 
                                        t]))]))/((2 - 
                                        y[1, t]) (1 + 
                                        Sqrt[(1 + 
                                        4 K P y[1, 
                                        t] (2 - y[1, t]))]))) 
                          1/P] /. x[3] -> x[3, t] /. x[4] -> x[4, t] /. 
                  x[1] -> x[1, t] /. x[2] -> x[2, t] /. 
              T -> T[t], 
Y[1, t] = 
            Chop[y[1, t] + y[3, t]], 
Y[2, t] = 
            Chop[y[2, t] + y[3, t]], 
X[1, t + 1] = 
            Chop[X[1, t] + h (X[1, t] - Y[1, t])], 
X[
              2, t + 1] = Chop[X[2, t] + h (X[2, t] - Y[2, t])], 
          Print[t], t++}] // Timing

tble2 = Table[{X[1, t], X[2, t]}, {t, 0, 499}];

plt2 = 
      ListPlot[tble2, PlotRange -> {{0, 1}, {0, 1}}, 
        PlotStyle -> RGBColor[1, 0, 0], PlotJoined ->  True, 
        Frame -> True]

Show[plt1, plt2]

t = 0; h = 0.5; X[1, 0] = 0.2; X[2, 0] = 0.1; x[1, -1] = 0.1; 
    x[2, -1] = 0.2; x[3, -1] = 0.1; x[4, -1] = 0.8; 
    y[1, -1] = 0.5;

While[
        t < 400, 
{
sol = 
            Chop[FindRoot[{X[1, t] ==  
                    x[1] + x[3], 
X[2, t] ==  
                    x[2] + x[3], 
                  Gam1 Gam2 x[1] x[2] - 
                      Gam3 Gam4 x[3] x[4]/Keq ==  
                    0, 
x[1] + x[2] + x[3] + 
                      x[4] ==  
                    1, 
Psat1 Gam1 x[
                            1]/((1 + 
                                  Sqrt[(1 + 4 K Psat1)])/(1 + 
                                  Sqrt[(1 + 
                                        4 K P y[
                                        1] (2 - y[1]))])) + 
                      Psat2 Gam2 x[
                            2]/((2 (1 - y[1] + 
                                      Sqrt[(1 + 
                                        4 K P y[
                                        1] (2 - 
                                        y[1]))]))/((2 - 
                                      y[1]) (1 + 
                                      Sqrt[(1 + 
                                        4 K P y[
                                        1] (2 - y[1]))]))) + 
                      Psat3 Gam3 x[
                            3]/((2 (1 - y[1] + 
                                      Sqrt[(1 + 
                                        4 K P y[
                                        1] (2 - 
                                        y[1]))]))/((2 - 
                                      y[1]) (1 + 
                                      Sqrt[(1 + 
                                        4 K P y[
                                        1] (2 - y[1]))]))) + 
                      Psat4 Gam4 x[
                            4]/((2 (1 - y[1] + 
                                      Sqrt[(1 + 
                                        4 K P y[
                                        1] (2 - 
                                        y[1]))]))/((2 - 
                                      y[1]) (1 + 
                                      Sqrt[(1 + 
                                        4 K P y[
                                        1] (2 - 
                                        y[1]))]))) ==  
                    P, 
y[1] ==  
                    Psat1 Gam1 x[
                          1]/((1 + Sqrt[(1 + 4 K Psat1)])/(1 + 
                                Sqrt[(1 + 
                                      4 K P y[
                                        1] (2 - y[1]))])) 1/
                        P}, {T, 350}, {x[3], x[3, t - 1]}, {x[4], 
                  x[4, t - 1]}, {x[1], x[1, t - 1]}, {x[2], 
                  x[2, t - 1]}, {y[1], y[1, t - 1]}, 
                MaxIterations -> 100000]], 
T[t] = 
            Chop[sol[[1, 2]]], 
x[3, t] = 
            Chop[sol[[2, 2]]], 
          x[4, t] = Chop[sol[[3, 2]]], 
x[1, t] = 
            Chop[sol[[4, 2]]], 
x[2, t] = 
            Chop[sol[[5, 2]]], 
y[1, t] = 
            Chop[sol[[6, 2]]], 
y[2, 
              t] = Chop[
                        Psat2 Gam2 x[2, 
                              t]/((2 (1 - y[1, t] + 
                                        Sqrt[(1 + 
                                        4 K P y[1, 
                                        t] (2 - y[1, 
                                        t]))]))/((2 - 
                                        y[1, t]) (1 + 
                                        Sqrt[(1 + 
                                        4 K P y[1, 
                                        t] (2 - y[1, t]))]))) 1/
                            P] /. x[3] -> x[3, t] /. x[4] -> x[4, t] /. 
                  x[1] -> x[1, t] /. x[2] -> x[2, t] /. 
              T -> T[t], 
y[3, 
              t] = Chop[
                        Psat3 Gam3 x[3, 
                              t]/((2 (1 - y[1, t] + 
                                        Sqrt[(1 + 
                                        4 K P y[1, 
                                        t] (2 - y[1, 
                                        t]))]))/((2 - 
                                        y[1, t]) (1 + 
                                        Sqrt[(1 + 
                                        4 K P y[1, 
                                        t] (2 - y[1, t]))]))) 
                          1/P] /. x[3] -> x[3, t] /. x[4] -> x[4, t] /. 
                  x[1] -> x[1, t] /. x[2] -> x[2, t] /. 
              T -> T[t], 
Y[1, t] = 
            Chop[y[1, t] + y[3, t]], 
Y[2, t] = 
            Chop[y[2, t] + y[3, t]], 
X[1, t + 1] = 
            Chop[X[1, t] - h (X[1, t] - Y[1, t])], 
X[
              2, t + 1] = Chop[X[2, t] - h (X[2, t] - Y[2, t])], 
          Print[t], t++}] // Timing

tble3 = Table[{X[1, t], X[2, t]}, {t, 0, 399}];

plt3 = 
      ListPlot[tble3, PlotRange -> {{0, 1}, {0, 1}}, 
        PlotStyle -> RGBColor[0, 0, 1], PlotJoined ->  True, 
        Frame -> True]

t = 0; h = 0.1; X[1, 0] = 0.2; X[2, 0] = 0.1; x[1, -1] = 0.1; 
    x[2, -1] = 0.2; x[3, -1] = 0.1; x[4, -1] = 0.8; 
    y[1, -1] = 0.5;

While[
        t < 200, 
{
sol = 
            Chop[FindRoot[{X[1, t] ==  
                    x[1] + x[3], 
X[2, t] ==  
                    x[2] + x[3], 
                  Gam1 Gam2 x[1] x[2] - 
                      Gam3 Gam4 x[3] x[4]/Keq ==  
                    0, 
x[1] + x[2] + x[3] + 
                      x[4] ==  
                    1, 
Psat1 Gam1 x[
                            1]/((1 + 
                                  Sqrt[(1 + 4 K Psat1)])/(1 + 
                                  Sqrt[(1 + 
                                        4 K P y[
                                        1] (2 - y[1]))])) + 
                      Psat2 Gam2 x[
                            2]/((2 (1 - y[1] + 
                                      Sqrt[(1 + 
                                        4 K P y[
                                        1] (2 - 
                                        y[1]))]))/((2 - 
                                      y[1]) (1 + 
                                      Sqrt[(1 + 
                                        4 K P y[
                                        1] (2 - y[1]))]))) + 
                      Psat3 Gam3 x[
                            3]/((2 (1 - y[1] + 
                                      Sqrt[(1 + 
                                        4 K P y[
                                        1] (2 - 
                                        y[1]))]))/((2 - 
                                      y[1]) (1 + 
                                      Sqrt[(1 + 
                                        4 K P y[
                                        1] (2 - y[1]))]))) + 
                      Psat4 Gam4 x[
                            4]/((2 (1 - y[1] + 
                                      Sqrt[(1 + 
                                        4 K P y[
                                        1] (2 - 
                                        y[1]))]))/((2 - 
                                      y[1]) (1 + 
                                      Sqrt[(1 + 
                                        4 K P y[
                                        1] (2 - 
                                        y[1]))]))) ==  
                    P, 
y[1] ==  
                    Psat1 Gam1 x[
                          1]/((1 + Sqrt[(1 + 4 K Psat1)])/(1 + 
                                Sqrt[(1 + 
                                      4 K P y[
                                        1] (2 - y[1]))])) 1/
                        P}, {T, 350}, {x[3], x[3, t - 1]}, {x[4], 
                  x[4, t - 1]}, {x[1], x[1, t - 1]}, {x[2], 
                  x[2, t - 1]}, {y[1], y[1, t - 1]}, 
                MaxIterations -> 100000]], 
T[t] = 
            Chop[sol[[1, 2]]], 
x[3, t] = 
            Chop[sol[[2, 2]]], 
          x[4, t] = Chop[sol[[3, 2]]], 
x[1, t] = 
            Chop[sol[[4, 2]]], 
x[2, t] = 
            Chop[sol[[5, 2]]], 
y[1, t] = 
            Chop[sol[[6, 2]]], 
y[2, 
              t] = Chop[
                        Psat2 Gam2 x[2, 
                              t]/((2 (1 - y[1, t] + 
                                        Sqrt[(1 + 
                                        4 K P y[1, 
                                        t] (2 - y[1, 
                                        t]))]))/((2 - 
                                        y[1, t]) (1 + 
                                        Sqrt[(1 + 
                                        4 K P y[1, 
                                        t] (2 - y[1, t]))]))) 1/
                            P] /. x[3] -> x[3, t] /. x[4] -> x[4, t] /. 
                  x[1] -> x[1, t] /. x[2] -> x[2, t] /. 
              T -> T[t], 
y[3, 
              t] = Chop[
                        Psat3 Gam3 x[3, 
                              t]/((2 (1 - y[1, t] + 
                                        Sqrt[(1 + 
                                        4 K P y[1, 
                                        t] (2 - y[1, 
                                        t]))]))/((2 - 
                                        y[1, t]) (1 + 
                                        Sqrt[(1 + 
                                        4 K P y[1, 
                                        t] (2 - y[1, t]))]))) 
                          1/P] /. x[3] -> x[3, t] /. x[4] -> x[4, t] /. 
                  x[1] -> x[1, t] /. x[2] -> x[2, t] /. 
              T -> T[t], 
Y[1, t] = 
            Chop[y[1, t] + y[3, t]], 
Y[2, t] = 
            Chop[y[2, t] + y[3, t]], 
X[1, t + 1] = 
            Chop[X[1, t] + h (X[1, t] - Y[1, t])], 
X[
              2, t + 1] = Chop[X[2, t] + h (X[2, t] - Y[2, t])], 
          Print[t], t++}] // Timing

tble4 = Table[{X[1, t], X[2, t]}, {t, 0, 199}];

plt4 = 
      ListPlot[tble4, PlotRange -> {{0, 1}, {0, 1}}, 
        PlotStyle -> RGBColor[0, 0, 1], PlotJoined ->  True, 
        Frame -> True]

Show[plt3, plt4]

t = 0; h = 0.5; X[1, 0] = 0.1; X[2, 0] = 0.95; x[1, -1] = 0.1; 
    x[2, -1] = 0.2; x[3, -1] = 0.1; x[4, -1] = 0.8; 
    y[1, -1] = 0.5;

While[
        t < 400, 
{
sol = 
            Chop[FindRoot[{X[1, t] ==  
                    x[1] + x[3], 
X[2, t] ==  
                    x[2] + x[3], 
                  Gam1 Gam2 x[1] x[2] - 
                      Gam3 Gam4 x[3] x[4]/Keq ==  
                    0, 
x[1] + x[2] + x[3] + 
                      x[4] ==  
                    1, 
Psat1 Gam1 x[
                            1]/((1 + 
                                  Sqrt[(1 + 4 K Psat1)])/(1 + 
                                  Sqrt[(1 + 
                                        4 K P y[
                                        1] (2 - y[1]))])) + 
                      Psat2 Gam2 x[
                            2]/((2 (1 - y[1] + 
                                      Sqrt[(1 + 
                                        4 K P y[
                                        1] (2 - 
                                        y[1]))]))/((2 - 
                                      y[1]) (1 + 
                                      Sqrt[(1 + 
                                        4 K P y[
                                        1] (2 - y[1]))]))) + 
                      Psat3 Gam3 x[
                            3]/((2 (1 - y[1] + 
                                      Sqrt[(1 + 
                                        4 K P y[
                                        1] (2 - 
                                        y[1]))]))/((2 - 
                                      y[1]) (1 + 
                                      Sqrt[(1 + 
                                        4 K P y[
                                        1] (2 - y[1]))]))) + 
                      Psat4 Gam4 x[
                            4]/((2 (1 - y[1] + 
                                      Sqrt[(1 + 
                                        4 K P y[
                                        1] (2 - 
                                        y[1]))]))/((2 - 
                                      y[1]) (1 + 
                                      Sqrt[(1 + 
                                        4 K P y[
                                        1] (2 - 
                                        y[1]))]))) ==  
                    P, 
y[1] ==  
                    Psat1 Gam1 x[
                          1]/((1 + Sqrt[(1 + 4 K Psat1)])/(1 + 
                                Sqrt[(1 + 
                                      4 K P y[
                                        1] (2 - y[1]))])) 1/
                        P}, {T, 350}, {x[3], x[3, t - 1]}, {x[4], 
                  x[4, t - 1]}, {x[1], x[1, t - 1]}, {x[2], 
                  x[2, t - 1]}, {y[1], y[1, t - 1]}, 
                MaxIterations -> 100000]], 
T[t] = 
            Chop[sol[[1, 2]]], 
x[3, t] = 
            Chop[sol[[2, 2]]], 
          x[4, t] = Chop[sol[[3, 2]]], 
x[1, t] = 
            Chop[sol[[4, 2]]], 
x[2, t] = 
            Chop[sol[[5, 2]]], 
y[1, t] = 
            Chop[sol[[6, 2]]], 
y[2, 
              t] = Chop[
                        Psat2 Gam2 x[2, 
                              t]/((2 (1 - y[1, t] + 
                                        Sqrt[(1 + 
                                        4 K P y[1, 
                                        t] (2 - y[1, 
                                        t]))]))/((2 - 
                                        y[1, t]) (1 + 
                                        Sqrt[(1 + 
                                        4 K P y[1, 
                                        t] (2 - y[1, t]))]))) 1/
                            P] /. x[3] -> x[3, t] /. x[4] -> x[4, t] /. 
                  x[1] -> x[1, t] /. x[2] -> x[2, t] /. 
              T -> T[t], 
y[3, 
              t] = Chop[
                        Psat3 Gam3 x[3, 
                              t]/((2 (1 - y[1, t] + 
                                        Sqrt[(1 + 
                                        4 K P y[1, 
                                        t] (2 - y[1, 
                                        t]))]))/((2 - 
                                        y[1, t]) (1 + 
                                        Sqrt[(1 + 
                                        4 K P y[1, 
                                        t] (2 - y[1, t]))]))) 
                          1/P] /. x[3] -> x[3, t] /. x[4] -> x[4, t] /. 
                  x[1] -> x[1, t] /. x[2] -> x[2, t] /. 
              T -> T[t], 
Y[1, t] = 
            Chop[y[1, t] + y[3, t]], 
Y[2, t] = 
            Chop[y[2, t] + y[3, t]], 
X[1, t + 1] = 
            Chop[X[1, t] - h (X[1, t] - Y[1, t])], 
X[
              2, t + 1] = Chop[X[2, t] - h (X[2, t] - Y[2, t])], 
          Print[t], t++}] // Timing

tble5 = Table[{X[1, t], X[2, t]}, {t, 0, 399}];

plt5 = 
      ListPlot[tble5, PlotRange -> {{0, 1}, {0, 1}}, 
        PlotStyle -> RGBColor[0, 1, 0], PlotJoined ->  True, 
        Frame -> True]

t = 0; h = 0.1; X[1, 0] = 0.1; X[2, 0] = 0.95; x[1, -1] = 0.1; 
    x[2, -1] = 0.2; x[3, -1] = 0.1; x[4, -1] = 0.8; 
    y[1, -1] = 0.5;

While[
        t < 200, 
{
sol = 
            Chop[FindRoot[{X[1, t] ==  
                    x[1] + x[3], 
X[2, t] ==  
                    x[2] + x[3], 
                  Gam1 Gam2 x[1] x[2] - 
                      Gam3 Gam4 x[3] x[4]/Keq ==  
                    0, 
x[1] + x[2] + x[3] + 
                      x[4] ==  
                    1, 
Psat1 Gam1 x[
                            1]/((1 + 
                                  Sqrt[(1 + 4 K Psat1)])/(1 + 
                                  Sqrt[(1 + 
                                        4 K P y[
                                        1] (2 - y[1]))])) + 
                      Psat2 Gam2 x[
                            2]/((2 (1 - y[1] + 
                                      Sqrt[(1 + 
                                        4 K P y[
                                        1] (2 - 
                                        y[1]))]))/((2 - 
                                      y[1]) (1 + 
                                      Sqrt[(1 + 
                                        4 K P y[
                                        1] (2 - y[1]))]))) + 
                      Psat3 Gam3 x[
                            3]/((2 (1 - y[1] + 
                                      Sqrt[(1 + 
                                        4 K P y[
                                        1] (2 - 
                                        y[1]))]))/((2 - 
                                      y[1]) (1 + 
                                      Sqrt[(1 + 
                                        4 K P y[
                                        1] (2 - y[1]))]))) + 
                      Psat4 Gam4 x[
                            4]/((2 (1 - y[1] + 
                                      Sqrt[(1 + 
                                        4 K P y[
                                        1] (2 - 
                                        y[1]))]))/((2 - 
                                      y[1]) (1 + 
                                      Sqrt[(1 + 
                                        4 K P y[
                                        1] (2 - 
                                        y[1]))]))) ==  
                    P, 
y[1] ==  
                    Psat1 Gam1 x[
                          1]/((1 + Sqrt[(1 + 4 K Psat1)])/(1 + 
                                Sqrt[(1 + 
                                      4 K P y[
                                        1] (2 - y[1]))])) 1/
                        P}, {T, 350}, {x[3], x[3, t - 1]}, {x[4], 
                  x[4, t - 1]}, {x[1], x[1, t - 1]}, {x[2], 
                  x[2, t - 1]}, {y[1], y[1, t - 1]}, 
                MaxIterations -> 100000]], 
T[t] = 
            Chop[sol[[1, 2]]], 
x[3, t] = 
            Chop[sol[[2, 2]]], 
          x[4, t] = Chop[sol[[3, 2]]], 
x[1, t] = 
            Chop[sol[[4, 2]]], 
x[2, t] = 
            Chop[sol[[5, 2]]], 
y[1, t] = 
            Chop[sol[[6, 2]]], 
y[2, 
              t] = Chop[
                        Psat2 Gam2 x[2, 
                              t]/((2 (1 - y[1, t] + 
                                        Sqrt[(1 + 
                                        4 K P y[1, 
                                        t] (2 - y[1, 
                                        t]))]))/((2 - 
                                        y[1, t]) (1 + 
                                        Sqrt[(1 + 
                                        4 K P y[1, 
                                        t] (2 - y[1, t]))]))) 1/
                            P] /. x[3] -> x[3, t] /. x[4] -> x[4, t] /. 
                  x[1] -> x[1, t] /. x[2] -> x[2, t] /. 
              T -> T[t], 
y[3, 
              t] = Chop[
                        Psat3 Gam3 x[3, 
                              t]/((2 (1 - y[1, t] + 
                                        Sqrt[(1 + 
                                        4 K P y[1, 
                                        t] (2 - y[1, 
                                        t]))]))/((2 - 
                                        y[1, t]) (1 + 
                                        Sqrt[(1 + 
                                        4 K P y[1, 
                                        t] (2 - y[1, t]))]))) 
                          1/P] /. x[3] -> x[3, t] /. x[4] -> x[4, t] /. 
                  x[1] -> x[1, t] /. x[2] -> x[2, t] /. 
              T -> T[t], 
Y[1, t] = 
            Chop[y[1, t] + y[3, t]], 
Y[2, t] = 
            Chop[y[2, t] + y[3, t]], 
X[1, t + 1] = 
            Chop[X[1, t] + h (X[1, t] - Y[1, t])], 
X[
              2, t + 1] = Chop[X[2, t] + h (X[2, t] - Y[2, t])], 
          Print[t], t++}] // Timing

tble6 = Table[{X[1, t], X[2, t]}, {t, 0, 199}];

plt6 = 
      ListPlot[tble6, PlotRange -> {{0, 1}, {0, 1}}, 
        PlotStyle -> RGBColor[0, 1, 0], PlotJoined ->  True, 
        Frame -> True]

Show[plt5, plt6]

Show[plt4, plt1, plt2, plt3, plt5, plt6, 
      Graphics[{Text[Butanol, {0, 1.1}], Text[Water, {-0.1, -0.1}], 
          Text[Acetic Acid, {1.1, -0.1}], 
          Text["Butyl Acetate", {1.1, 1.1}]}], AspectRatio -> 1, 
      PlotRange -> All, Frame -> False, 
      Epilog -> {Line[{{0, 1}, {1, 1}}], 
          Line[{{1, 0}, {1, 1}}]}]

The coordinates of the reactive azeotrope X1=0.04238 and X2=0.218399.

Last[tble1]

Last[tble3]

Last[tble5]