Residue Curve Map for Methyl Acetate Chemistry 
Author's Data: Housam BINOUS
Department of Chemical Engineering
National Institute of Applied Sciences and Technology
Tunis, TUNISIA
Email: binoushousam@yahoo.com

Activity Coefficient using Wilson model
Component 1: Acetic Acid
Component 2: Methanol
Component 3: Methyl Acetate
Component 4: Water

G[i_] := Exp[
                1 - Log[(∑_(j = 1)^4) x[j] A[i, 
                          j]] - (∑_(k = 1)^4)(((x[k] A[k, i])/((∑_(j = 1)^4) x[j] A[k, j])))] /. x[3] -> x3 /. 
            x[4] -> x4 /. x[1] -> x1 /. 
        x[2] -> x2

A[i_, j_] := V[j]/V[i] Exp[-a[i, j]/(R T)]

a[1, 1] = 0;
a[1, 2] = 2535.2019;
a[1, 3] = 1123.1444;
a[1, 4] = 237.5248;
a[2, 1] = -547.5248;
a[2, 2] = 0;
a[2, 3] = 813.1843;
a[2, 4] = 107.3832;
a[3, 1] = -696.5031;
a[3, 2] = -31.1932;
a[3, 3] = 0;
a[3, 4] = 645.7225;
a[4, 1] = 658.0266;
a[4, 2] = 469.5509;
a[4, 3] = 1918.232;
a[4, 4] = 0;

V[1] = 57.54;
V[2] = 44.44;
V[3] = 79.84;
V[4] = 18.07;

Gas Constant and Total Pressure

R = 1.987; P = 1.013 10^5;

Vapor Pressure Using Antoine Equation and Normal Boiling Temperatures

Psat[i_] := Exp[(Aa[i] + Ba[i]/(T + Ca[i]))]

Aa[1] = 22.1001; Aa[2] = 23.49989;
Aa[3] = 21.1520; Aa[4] = 23.2256;
Ba[1] = -3654.62; Ba[2] = -3643.31362; Ba[3] = -2662.78; 
    Ba[4] = -3835.18;
Ca[1] = -45.392; Ca[2] = -33.434; Ca[3] = -53.460; 
    Ca[4] = -45.343;

FindRoot[Psat[1] ==  1.013 10^5, {T, 300}]

FindRoot[Psat[2] ==  1.013 10^5, {T, 300}]

FindRoot[Psat[3] ==  1.013 10^5, {T, 300}]

FindRoot[Psat[4] ==  1.013 10^5, {T, 300}]

Equilibrium Constant for the esterification of acetic acid with methanol to produce methyl acetate and water:
	C(Subscript[H, 3])COOH + C(Subscript[H, 3])OH ⇌ C(Subscript[H, 3])COOC(Subscript[H, 3]) +  (Subscript[H, 2])O

Keq = Exp[0.83983 + 782.98/T];

Equilibrium Constant for the Dimerization of Acetic Acid

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

Computing the Residue Curve Map takes about 15 minutes on a Pentium 800 MHz PC. Azeotrope compositions is : X1=0.669, X2=1 ( transformed compositions ).

t = 0; h = 0.1; X[1, 0] = 0.2; X[2, 0] = 0.8; 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 < 232, 
{
sol = 
            Chop[FindRoot[{X[1, t] ==  
                    x1 + x3, 
X[2, t] ==  x2 + x3, 
                  G[1] G[2] x1 x2 - G[3] G[4] x3 x4/Keq ==  
                    0, 
x1 + x2 + x3 + x4 ==  
                    1, 
Psat[1] G[1] 
                        x1/((1 + Sqrt[(1 + 4 K Psat[1])])/(1 + 
                                  Sqrt[(1 + 
                                        4 K P y1 (2 - y1))])) + 
                      Psat[2] G[2] 
                        x2/((2 (1 - y1 + 
                                      Sqrt[(1 + 
                                        4 K P y1 (2 - 
                                        y1))]))/((2 - 
                                      y1) (1 + 
                                      Sqrt[(1 + 
                                        4 K P y1 (2 - 
                                        y1))]))) + 
                      Psat[3] G[3] 
                        x3/((2 (1 - y1 + 
                                      Sqrt[(1 + 
                                        4 K P y1 (2 - 
                                        y1))]))/((2 - 
                                      y1) (1 + 
                                      Sqrt[(1 + 
                                        4 K P y1 (2 - 
                                        y1))]))) + 
                      Psat[4] G[4] 
                        x4/((2 (1 - y1 + 
                                      Sqrt[(1 + 
                                        4 K P y1 (2 - 
                                        y1))]))/((2 - 
                                      y1) (1 + 
                                      Sqrt[(1 + 
                                        4 K P y1 (2 - 
                                        y1))]))) ==  
                    P, 
y1 ==  
                    Psat[1] G[1] 
                      x1/((1 + Sqrt[(1 + 4 K Psat[1])])/(1 + 
                                Sqrt[(1 + 
                                      4 K P y1 (2 - y1))])) 1/
                        P}, {T, 350}, {x3, x[3, t - 1]}, {x4, 
                  x[4, t - 1]}, {x1, x[1, t - 1]}, {x2, x[2, t - 1]}, {y1, 
                  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[
                        Psat[2] G[2] 
                          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] /. x3 -> x[3, t] /. x4 -> x[4, t] /. 
                  x1 -> x[1, t] /. x2 -> x[2, t] /. 
              T -> T[t], 
y[3, 
              t] = Chop[
                        Psat[3] G[3] 
                          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] /. x3 -> x[3, t] /. x4 -> x[4, t] /. 
                  x1 -> x[1, t] /. x2 -> 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])], t++}] // 
      Timing

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

Last[tble1]

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

t = 0; h = 0.01; X[1, 0] = 0.2; X[2, 0] = 0.8; 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] ==  
                    x1 + x3, 
X[2, t] ==  x2 + x3, 
                  G[1] G[2] x1 x2 - G[3] G[4] x3 x4/Keq ==  
                    0, 
x1 + x2 + x3 + x4 ==  
                    1, 
Psat[1] G[1] 
                        x1/((1 + Sqrt[(1 + 4 K Psat[1])])/(1 + 
                                  Sqrt[(1 + 
                                        4 K P y1 (2 - y1))])) + 
                      Psat[2] G[2] 
                        x2/((2 (1 - y1 + 
                                      Sqrt[(1 + 
                                        4 K P y1 (2 - 
                                        y1))]))/((2 - 
                                      y1) (1 + 
                                      Sqrt[(1 + 
                                        4 K P y1 (2 - 
                                        y1))]))) + 
                      Psat[3] G[3] 
                        x3/((2 (1 - y1 + 
                                      Sqrt[(1 + 
                                        4 K P y1 (2 - 
                                        y1))]))/((2 - 
                                      y1) (1 + 
                                      Sqrt[(1 + 
                                        4 K P y1 (2 - 
                                        y1))]))) + 
                      Psat[4] G[4] 
                        x4/((2 (1 - y1 + 
                                      Sqrt[(1 + 
                                        4 K P y1 (2 - 
                                        y1))]))/((2 - 
                                      y1) (1 + 
                                      Sqrt[(1 + 
                                        4 K P y1 (2 - 
                                        y1))]))) ==  
                    P, 
y1 ==  
                    Psat[1] G[1] 
                      x1/((1 + Sqrt[(1 + 4 K Psat[1])])/(1 + 
                                Sqrt[(1 + 
                                      4 K P y1 (2 - y1))])) 1/
                        P}, {T, 350}, {x3, x[3, t - 1]}, {x4, 
                  x[4, t - 1]}, {x1, x[1, t - 1]}, {x2, x[2, t - 1]}, {y1, 
                  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[
                        Psat[2] G[2] 
                          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] /. x3 -> x[3, t] /. x4 -> x[4, t] /. 
                  x1 -> x[1, t] /. x2 -> x[2, t] /. 
              T -> T[t], 
y[3, 
              t] = Chop[
                        Psat[3] G[3] 
                          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] /. x3 -> x[3, t] /. x4 -> x[4, t] /. 
                  x1 -> x[1, t] /. x2 -> 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])], t++}] // 
      Timing

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

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

t = 0; h = 0.1; X[1, 0] = 0.5; X[2, 0] = 0.5;

While[
        t < 242, 
{
sol = 
            Chop[FindRoot[{X[1, t] ==  
                    x1 + x3, 
X[2, t] ==  x2 + x3, 
                  G[1] G[2] x1 x2 - G[3] G[4] x3 x4/Keq ==  
                    0 ==  
                    0, 
x1 + x2 + x3 + x4 ==  
                    1, 
Psat[1] G[1] 
                        x1/((1 + Sqrt[(1 + 4 K Psat[1])])/(1 + 
                                  Sqrt[(1 + 
                                        4 K P y1 (2 - y1))])) + 
                      Psat[2] G[2] 
                        x2/((2 (1 - y1 + 
                                      Sqrt[(1 + 
                                        4 K P y1 (2 - 
                                        y1))]))/((2 - 
                                      y1) (1 + 
                                      Sqrt[(1 + 
                                        4 K P y1 (2 - 
                                        y1))]))) + 
                      Psat[3] G[3] 
                        x3/((2 (1 - y1 + 
                                      Sqrt[(1 + 
                                        4 K P y1 (2 - 
                                        y1))]))/((2 - 
                                      y1) (1 + 
                                      Sqrt[(1 + 
                                        4 K P y1 (2 - 
                                        y1))]))) + 
                      Psat[4] G[4] 
                        x4/((2 (1 - y1 + 
                                      Sqrt[(1 + 
                                        4 K P y1 (2 - 
                                        y1))]))/((2 - 
                                      y1) (1 + 
                                      Sqrt[(1 + 
                                        4 K P y1 (2 - 
                                        y1))]))) ==  
                    P, 
y1 ==  
                    Psat[1] G[1] 
                      x1/((1 + Sqrt[(1 + 4 K Psat[1])])/(1 + 
                                Sqrt[(1 + 
                                      4 K P y1 (2 - y1))])) 1/
                        P}, {T, 350}, {x3, x[3, t - 1]}, {x4, 
                  x[4, t - 1]}, {x1, x[1, t - 1]}, {x2, x[2, t - 1]}, {y1, 
                  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[
                        Psat[2] G[2] 
                          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] /. x3 -> x[3, t] /. x4 -> x[4, t] /. 
                  x1 -> x[1, t] /. x2 -> x[2, t] /. 
              T -> T[t], 
y[3, 
              t] = Chop[
                        Psat[3] G[3] 
                          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] /. x3 -> x[3, t] /. x4 -> x[4, t] /. 
                  x1 -> x[1, t] /. x2 -> 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])], t++}] // 
      Timing

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

Last[tble3]

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

t = 0; h = 0.05; X[1, 0] = 0.5; X[2, 0] = 0.5; 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 < 320, 
{
sol = 
            Chop[FindRoot[{X[1, t] ==  
                    x1 + x3, 
X[2, t] ==  x2 + x3, 
                  G[1] G[2] x1 x2 - G[3] G[4] x3 x4/Keq ==  
                    0, 
x1 + x2 + x3 + x4 ==  
                    1, 
Psat[1] G[1] 
                        x1/((1 + Sqrt[(1 + 4 K Psat[1])])/(1 + 
                                  Sqrt[(1 + 
                                        4 K P y1 (2 - y1))])) + 
                      Psat[2] G[2] 
                        x2/((2 (1 - y1 + 
                                      Sqrt[(1 + 
                                        4 K P y1 (2 - 
                                        y1))]))/((2 - 
                                      y1) (1 + 
                                      Sqrt[(1 + 
                                        4 K P y1 (2 - 
                                        y1))]))) + 
                      Psat[3] G[3] 
                        x3/((2 (1 - y1 + 
                                      Sqrt[(1 + 
                                        4 K P y1 (2 - 
                                        y1))]))/((2 - 
                                      y1) (1 + 
                                      Sqrt[(1 + 
                                        4 K P y1 (2 - 
                                        y1))]))) + 
                      Psat[4] G[4] 
                        x4/((2 (1 - y1 + 
                                      Sqrt[(1 + 
                                        4 K P y1 (2 - 
                                        y1))]))/((2 - 
                                      y1) (1 + 
                                      Sqrt[(1 + 
                                        4 K P y1 (2 - 
                                        y1))]))) ==  
                    P, 
y1 ==  
                    Psat[1] G[1] 
                      x1/((1 + Sqrt[(1 + 4 K Psat[1])])/(1 + 
                                Sqrt[(1 + 
                                      4 K P y1 (2 - y1))])) 1/
                        P}, {T, 350}, {x3, x[3, t - 1]}, {x4, 
                  x[4, t - 1]}, {x1, x[1, t - 1]}, {x2, x[2, t - 1]}, {y1, 
                  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[
                        Psat[2] G[2] 
                          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] /. x3 -> x[3, t] /. x4 -> x[4, t] /. 
                  x1 -> x[1, t] /. x2 -> x[2, t] /. 
              T -> T[t], 
y[3, 
              t] = Chop[
                        Psat[3] G[3] 
                          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] /. x3 -> x[3, t] /. x4 -> x[4, t] /. 
                  x1 -> x[1, t] /. x2 -> 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])], t++}] // 
      Timing

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

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

t = 0; h = 0.1; X[1, 0] = 0.8; X[2, 0] = 0.8;

While[
        t < 228, 
{
sol = 
            Chop[FindRoot[{X[1, t] ==  
                    x1 + x3, 
X[2, t] ==  x2 + x3, 
                  G[1] G[2] x1 x2 - G[3] G[4] x3 x4/Keq ==  
                    0, 
x1 + x2 + x3 + x4 ==  
                    1, 
Psat[1] G[1] 
                        x1/((1 + Sqrt[(1 + 4 K Psat[1])])/(1 + 
                                  Sqrt[(1 + 
                                        4 K P y1 (2 - y1))])) + 
                      Psat[2] G[2] 
                        x2/((2 (1 - y1 + 
                                      Sqrt[(1 + 
                                        4 K P y1 (2 - 
                                        y1))]))/((2 - 
                                      y1) (1 + 
                                      Sqrt[(1 + 
                                        4 K P y1 (2 - 
                                        y1))]))) + 
                      Psat[3] G[3] 
                        x3/((2 (1 - y1 + 
                                      Sqrt[(1 + 
                                        4 K P y1 (2 - 
                                        y1))]))/((2 - 
                                      y1) (1 + 
                                      Sqrt[(1 + 
                                        4 K P y1 (2 - 
                                        y1))]))) + 
                      Psat[4] G[4] 
                        x4/((2 (1 - y1 + 
                                      Sqrt[(1 + 
                                        4 K P y1 (2 - 
                                        y1))]))/((2 - 
                                      y1) (1 + 
                                      Sqrt[(1 + 
                                        4 K P y1 (2 - 
                                        y1))]))) ==  
                    P, 
y1 ==  
                    Psat[1] G[1] 
                      x1/((1 + Sqrt[(1 + 4 K Psat[1])])/(1 + 
                                Sqrt[(1 + 
                                      4 K P y1 (2 - y1))])) 1/
                        P}, {T, 350}, {x3, x[3, t - 1]}, {x4, 
                  x[4, t - 1]}, {x1, x[1, t - 1]}, {x2, x[2, t - 1]}, {y1, 
                  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[
                        Psat[2] G[2] 
                          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] /. x3 -> x[3, t] /. x4 -> x[4, t] /. 
                  x1 -> x[1, t] /. x2 -> x[2, t] /. 
              T -> T[t], 
y[3, 
              t] = Chop[
                        Psat[3] G[3] 
                          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] /. x3 -> x[3, t] /. x4 -> x[4, t] /. 
                  x1 -> x[1, t] /. x2 -> 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])], t++}] // 
      Timing

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

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

t = 0; h = 0.05; X[1, 0] = 0.8; X[2, 0] = 0.8; 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 < 276, 
{
sol = 
            Chop[FindRoot[{X[1, t] ==  
                    x1 + x3, 
X[2, t] ==  x2 + x3, 
                  G[1] G[2] x1 x2 - G[3] G[4] x3 x4/Keq ==  
                    0, 
x1 + x2 + x3 + x4 ==  
                    1, 
Psat[1] G[1] 
                        x1/((1 + Sqrt[(1 + 4 K Psat[1])])/(1 + 
                                  Sqrt[(1 + 
                                        4 K P y1 (2 - y1))])) + 
                      Psat[2] G[2] 
                        x2/((2 (1 - y1 + 
                                      Sqrt[(1 + 
                                        4 K P y1 (2 - 
                                        y1))]))/((2 - 
                                      y1) (1 + 
                                      Sqrt[(1 + 
                                        4 K P y1 (2 - 
                                        y1))]))) + 
                      Psat[3] G[3] 
                        x3/((2 (1 - y1 + 
                                      Sqrt[(1 + 
                                        4 K P y1 (2 - 
                                        y1))]))/((2 - 
                                      y1) (1 + 
                                      Sqrt[(1 + 
                                        4 K P y1 (2 - 
                                        y1))]))) + 
                      Psat[4] G[4] 
                        x4/((2 (1 - y1 + 
                                      Sqrt[(1 + 
                                        4 K P y1 (2 - 
                                        y1))]))/((2 - 
                                      y1) (1 + 
                                      Sqrt[(1 + 
                                        4 K P y1 (2 - 
                                        y1))]))) ==  
                    P, 
y1 ==  
                    Psat[1] G[1] 
                      x1/((1 + Sqrt[(1 + 4 K Psat[1])])/(1 + 
                                Sqrt[(1 + 
                                      4 K P y1 (2 - y1))])) 1/
                        P}, {T, 350}, {x3, x[3, t - 1]}, {x4, 
                  x[4, t - 1]}, {x1, x[1, t - 1]}, {x2, x[2, t - 1]}, {y1, 
                  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[
                        Psat[2] G[2] 
                          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] /. x3 -> x[3, t] /. x4 -> x[4, t] /. 
                  x1 -> x[1, t] /. x2 -> x[2, t] /. 
              T -> T[t], 
y[3, 
              t] = Chop[
                        Psat[3] G[3] 
                          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] /. x3 -> x[3, t] /. x4 -> x[4, t] /. 
                  x1 -> x[1, t] /. x2 -> 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])], t++}] // 
      Timing

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

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

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