Residue Curve Map for Isopropyl 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 Coefficients based on NRTL model
Component 1: Acetic Acid
Component 2: Isopropanol
Component 3: Isopropyl Acetate
Component 4: Water

Gam[i_]:=Exp[(((∑_(j = 1)^4)( tau[j, i] G[j, i] x[j]))/((∑_(l = 1)^4) (G[l, i] x[l])))+(∑_(j = 1)^4) (((G[i, j] x[j])/((∑_(l = 1)^4)( G[l, j] x[l])))(tau[i, j]-(((∑_(r = 1)^4)( tau[r, j] 
                                    G[r, j] x[r]))/((∑_(l = 1)^4) (G[l, j] x[l])))))]/.x[1] -> x1/.x [2] -> x2/.x[3] -> x3/.x[4] -> x4

G[i_, j_] := Exp[-a[i, j] tau[i, j]]

tau[i_, j_] := g[i, j]/(R  T)

g[1, 1] = 0; g[2, 1] = 81.3926; g[3, 1] = 154.7885; 
    g[4, 1] = 842.6081;
g[1, 2] = -281.4482; g[2, 2] = 0; g[3, 2] = 140.0972; 
    g[4, 2] = 1655.2550;
g[1, 3] = 141.0082; g[2, 3] = 269.9609; g[3, 3] = 0; 
    g[4, 3] = 1270.2036;
g[1, 4] = -219.7238; g[2, 4] = 39.8541; g[3, 4] = 1165.709; 
    g[4, 4] = 0;

a[1, 1] = 0; a[2, 1] = 0.3048; a[3, 1] = 0.3014; 
    a[4, 1] = 0.2997;
a[1, 2] = 0.3048; a[2, 2] = 0; a[3, 2] = 0.3009; 
    a[4, 2] = 0.3255;
a[1, 3] = 0.3014; a[2, 3] = 0.3009; a[3, 3] = 0; 
    a[4, 3] = 0.3300;
a[1, 4] = 0.2997; a[2, 4] = 0.3255; a[3, 4] = 0.3300; 
    a[4, 4] = 0;

Gas Constant and Total Pressure

R = 1.987; P = 1.013 10^5;

Vapor Pressure based on Antoine Equation and Normal Boiling Points

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

Aa[1] = 23.3618; Aa[2] = 25.3358;
Aa[3] = 21.7798; Aa[4] = 23.4776;
Ba[1] = -4457.83; Ba[2] = -4628.96; Ba[3] = -3307.73; 
    Ba[4] = -3984.92;
Ca[1] = -14.699; Ca[2] = -20.514; Ca[3] = -39.485; 
    Ca[4] = -39.724;

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 Dimerization of Acetic Acid

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

Equilibrium Constant for the esterification of acetic acid with isopropanol to give isopropyl acetate and water:
	C(Subscript[H, 3])COOH + C(Subscript[H, 3])CH(OH)C(Subscript[H, 3]) ⇌ C(Subscript[H, 3])COOCH(C(Subscript[H, 3]))_2 + (Subscript[H, 2])O

Keq = 8.7;

Computation of Residue Curve Map takes approximately 5 minutes on a Pentium 800 MHz PC. Reactive Azeotrope compositions is : X1=0.231, X2=0.748 ( transformed compositions ) or x1=0.048, x2=0.565, x3=0.183, x4=0.204.

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 < 179, 
{
sol = 
            Chop[FindRoot[{X[1, t] ==  
                    x1 + x3, 
X[2, t] ==  x2 + x3, 
                  Gam[1] Gam[2] x1 x2 - 
                      Gam[3] Gam[4] x3 x4/Keq ==  
                    0, 
x1 + x2 + x3 + x4 ==  
                    1, 
Psat[1] Gam[1] 
                        x1/((1 + Sqrt[(1 + 4 K Psat[1])])/(1 + 
                                  Sqrt[(1 + 
                                        4 K P y1 (2 - y1))])) + 
                      Psat[2] Gam[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] Gam[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] Gam[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] Gam[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] Gam[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] Gam[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

x[1, 177]

x[2, 177]

x[3, 177]

x[4, 177]

y[1, 177]

y[2, 177]

y[3, 177]

{X[1,177],X[2,177]}

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

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

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 < 25, 
{
sol = 
            Chop[FindRoot[{X[1, t] ==  
                    x1 + x3, 
X[2, t] ==  x2 + x3, 
                  Gam[1] Gam[2] x1 x2 - 
                      Gam[3] Gam[4] x3 x4/Keq ==  
                    0, 
x1 + x2 + x3 + x4 ==  
                    1, 
Psat[1] Gam[1] 
                        x1/((1 + Sqrt[(1 + 4 K Psat[1])])/(1 + 
                                  Sqrt[(1 + 
                                        4 K P y1 (2 - y1))])) + 
                      Psat[2] Gam[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] Gam[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] Gam[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] Gam[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] Gam[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] Gam[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, 24}];

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

t = 0; h = 0.5; 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 < 28, 
{
sol = 
            Chop[FindRoot[{X[1, t] ==  
                    x1 + x3, 
X[2, t] ==  x2 + x3, 
                  Gam[1] Gam[2] x1 x2 - 
                      Gam[3] Gam[4] x3 x4/Keq ==  
                    0, 
x1 + x2 + x3 + x4 ==  
                    1, 
Psat[1] Gam[1] 
                        x1/((1 + Sqrt[(1 + 4 K Psat[1])])/(1 + 
                                  Sqrt[(1 + 
                                        4 K P y1 (2 - y1))])) + 
                      Psat[2] Gam[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] Gam[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] Gam[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] Gam[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] Gam[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] Gam[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, 27}];

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

t = 0; h = 0.5; 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 < 80, 
{
sol = 
            Chop[FindRoot[{X[1, t] ==  
                    x1 + x3, 
X[2, t] ==  x2 + x3, 
                  Gam[1] Gam[2] x1 x2 - 
                      Gam[3] Gam[4] x3 x4/Keq ==  
                    0, 
x1 + x2 + x3 + x4 ==  
                    1, 
Psat[1] Gam[1] 
                        x1/((1 + Sqrt[(1 + 4 K Psat[1])])/(1 + 
                                  Sqrt[(1 + 
                                        4 K P y1 (2 - y1))])) + 
                      Psat[2] Gam[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] Gam[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] Gam[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] Gam[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] Gam[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] Gam[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, 79}];

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

t = 0; h = 0.5; X[1, 0] = 0.2; X[2, 0] = 0.9; 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 < 51, 
{
sol = 
            Chop[FindRoot[{X[1, t] ==  
                    x1 + x3, 
X[2, t] ==  x2 + x3, 
                  Gam[1] Gam[2] x1 x2 - 
                      Gam[3] Gam[4] x3 x4/Keq ==  
                    0, 
x1 + x2 + x3 + x4 ==  
                    1, 
Psat[1] Gam[1] 
                        x1/((1 + Sqrt[(1 + 4 K Psat[1])])/(1 + 
                                  Sqrt[(1 + 
                                        4 K P y1 (2 - y1))])) + 
                      Psat[2] Gam[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] Gam[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] Gam[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] Gam[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] Gam[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] Gam[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, 50}];

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

t = 0; h = 0.5; X[1, 0] = 0.2; X[2, 0] = 0.9; 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 < 61, 
{
sol = 
            Chop[FindRoot[{X[1, t] ==  
                    x1 + x3, 
X[2, t] ==  x2 + x3, 
                  Gam[1] Gam[2] x1 x2 - 
                      Gam[3] Gam[4] x3 x4/Keq ==  
                    0, 
x1 + x2 + x3 + x4 ==  
                    1, 
Psat[1] Gam[1] 
                        x1/((1 + Sqrt[(1 + 4 K Psat[1])])/(1 + 
                                  Sqrt[(1 + 
                                        4 K P y1 (2 - y1))])) + 
                      Psat[2] Gam[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] Gam[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] Gam[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] Gam[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] Gam[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] Gam[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, 79}];

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

t = 0; h = 0.05; X[1, 0] = 0.1; X[2, 0] = 0.2; 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 < 300, 
{
sol = 
            Chop[FindRoot[{X[1, t] ==  
                    x1 + x3, 
X[2, t] ==  x2 + x3, 
                  Gam[1] Gam[2] x1 x2 - 
                      Gam[3] Gam[4] x3 x4/Keq ==  
                    0, 
x1 + x2 + x3 + x4 ==  
                    1, 
Psat[1] Gam[1] 
                        x1/((1 + Sqrt[(1 + 4 K Psat[1])])/(1 + 
                                  Sqrt[(1 + 
                                        4 K P y1 (2 - y1))])) + 
                      Psat[2] Gam[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] Gam[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] Gam[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] Gam[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] Gam[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] Gam[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

tble7 = Table[{X[1, t], X[2, t]}, {t, 0, 299}];

plt7 = 
      ListPlot[tble7, PlotRange -> {{0, 1}, {-0.1, 1}}, 
        PlotStyle -> RGBColor[1, 0, 1], PlotJoined ->  True, 
        Frame -> True]

t = 0; h = 0.5; X[1, 0] = 0.1; X[2, 0] = 0.2; 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 < 61, 
{
sol = 
            Chop[FindRoot[{X[1, t] ==  
                    x1 + x3, 
X[2, t] ==  x2 + x3, 
                  Gam[1] Gam[2] x1 x2 - 
                      Gam[3] Gam[4] x3 x4/Keq ==  
                    0, 
x1 + x2 + x3 + x4 ==  
                    1, 
Psat[1] Gam[1] 
                        x1/((1 + Sqrt[(1 + 4 K Psat[1])])/(1 + 
                                  Sqrt[(1 + 
                                        4 K P y1 (2 - y1))])) + 
                      Psat[2] Gam[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] Gam[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] Gam[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] Gam[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] Gam[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] Gam[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

tble8 = Table[{X[1, t], X[2, t]}, {t, 0, 60}];

plt8 = 
      ListPlot[tble8, PlotRange -> {{0, 1}, {-0.1, 1}}, 
        PlotStyle -> RGBColor[1, 0, 1], PlotJoined ->  True, 
        Frame -> True]

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