Simple distillation Residue Curve Map for Ternary Systems (A, B and C) 
with Reaction A+B=C
Author's Data: Housam BINOUS
Department of Chemical Engineering
National Institute of Applied Sciences and Technology
Tunis, TUNISIA
Email: binoushousam@yahoo.com

Total Pressure

P = 1.013 10^5;

Vapor Pressure Using Antoine Equation of the three components A,B,C of a hypothetical ideal mixture.

A1 = 23.1962; B1 = -3502.77; C1 = 0.0;
A2 = 22.6867; B2 = -3502.77; C2 = 0.0;
A3 = 21.5872; B3 = -3502.77; C3 = 0.0;

PS1[i_] := Exp[A1 + B1/(C1 + T)];
PS2[i_] := Exp[A2 + B2/(C2 + T)];
PS3[i_] := Exp[A3 + B3/(C3 + T)];

Normal Boiling Points of components A,B and C.

FindRoot[PS1[i] == 1.013 10^5, {T, 273}]

{T -> 300.142}

FindRoot[PS2[i] == 1.013 10^5, {T, 273}]

{T -> 313.844}

FindRoot[PS3[i] == 1.013 10^5, {T, 273}]

{T -> 348.141}

Residue Curve Map When Damkoler number (Da) equals 0
(i.e. no reaction case).

h = 0.1; x1[0] = 0.2; x2[0] = 0.2; x3[0] = 0.6; i = 0;

While[
        i < 200, {T = FindRoot[
                P == 
                  Evaluate[
                    x1[i] PS1[i] + x2[i] PS2[i] + x3[i] PS3[i] ], {T, 
                  273}][[1, 2]], 
          y1[i] = x1[i] PS1[i] /P, 
y2[i] = 
            x2[i] PS2[i] /P, 
y3[i] = x3[i] PS3[i] /P, 
          x1[i + 1] = x1[i] - h (x1[i] - y1[i]), 
          x2[i + 1] = x2[i] - h (x2[i] - y2[i]), 
          x3[i + 1] = x3[i] - h (x3[i] - y3[i]), i++}];

t1 = Table[{x1[i], x2[i]}, {i, 1, 200}];

plt1 = ListPlot[t1, PlotRange -> {{0, 1}, {0, 1}}];

h = 0.1; x1[0] = 0.2; x2[0] = 0.2; x3[0] = 0.6; i = 0;

While[
        i < 200, {T = FindRoot[
                P == 
                  Evaluate[
                    x1[i] PS1[i] + x2[i] PS2[i] + x3[i] PS3[i] ], {T, 
                  273}][[1, 2]], 
          y1[i] = x1[i] PS1[i] /P, 
y2[i] = 
            x2[i] PS2[i] /P, 
y3[i] = x3[i] PS3[i] /P, 
          x1[i + 1] = x1[i] + h (x1[i] - y1[i]), 
          x2[i + 1] = x2[i] + h (x2[i] - y2[i]), 
          x3[i + 1] = x3[i] + h (x3[i] - y3[i]), i++}];

t2 = Table[{x1[i], x2[i]}, {i, 1, 200}];

plt2 = ListPlot[t2, PlotRange -> {{0, 1}, {0, 1}}];

h = 0.1; x1[0] = 0.4; x2[0] = 0.4; x3[0] = 0.2; i = 0;

While[
        i < 200, {T = FindRoot[
                P == 
                  Evaluate[
                    x1[i] PS1[i] + x2[i] PS2[i] + x3[i] PS3[i] ], {T, 
                  273}][[1, 2]], 
          y1[i] = x1[i] PS1[i] /P, 
y2[i] = 
            x2[i] PS2[i] /P, 
y3[i] = x3[i] PS3[i] /P, 
          x1[i + 1] = x1[i] - h (x1[i] - y1[i]), 
          x2[i + 1] = x2[i] - h (x2[i] - y2[i]), 
          x3[i + 1] = x3[i] - h (x3[i] - y3[i]), i++}];

t3 = Table[{x1[i], x2[i]}, {i, 1, 200}];

plt3 = ListPlot[t3, PlotRange -> {{0, 1}, {0, 1}}];

h = 0.1; x1[0] = 0.4; x2[0] = 0.4; x3[0] = 0.2; i = 0;

While[
        i < 200, {T = FindRoot[
                P == 
                  Evaluate[
                    x1[i] PS1[i] + x2[i] PS2[i] + x3[i] PS3[i] ], {T, 
                  273}][[1, 2]], 
          y1[i] = x1[i] PS1[i] /P, 
y2[i] = 
            x2[i] PS2[i] /P, 
y3[i] = x3[i] PS3[i] /P, 
          x1[i + 1] = x1[i] + h (x1[i] - y1[i]), 
          x2[i + 1] = x2[i] + h (x2[i] - y2[i]), 
          x3[i + 1] = x3[i] + h (x3[i] - y3[i]), i++}];

t4 = Table[{x1[i], x2[i]}, {i, 1, 200}];

plt4 = ListPlot[t4, PlotRange -> {{0, 1}, {0, 1}}];

h = 0.1; x1[0] = 0.8; x2[0] = 0.1; x3[0] = 0.1; i = 0;

While[
        i < 200, {T = FindRoot[
                P == 
                  Evaluate[
                    x1[i] PS1[i] + x2[i] PS2[i] + x3[i] PS3[i] ], {T, 
                  273}][[1, 2]], 
          y1[i] = x1[i] PS1[i] /P, 
y2[i] = 
            x2[i] PS2[i] /P, 
y3[i] = x3[i] PS3[i] /P, 
          x1[i + 1] = x1[i] - h (x1[i] - y1[i]), 
          x2[i + 1] = x2[i] - h (x2[i] - y2[i]), 
          x3[i + 1] = x3[i] - h (x3[i] - y3[i]), i++}];

t5 = Table[{x1[i], x2[i]}, {i, 1, 200}];

plt5 = ListPlot[t5, PlotRange -> {{0, 1}, {0, 1}}];

h = 0.1; x1[0] = 0.8; x2[0] = 0.1; x3[0] = 0.1; i = 0;

While[
        i < 200, {T = FindRoot[
                P == 
                  Evaluate[
                    x1[i] PS1[i] + x2[i] PS2[i] + x3[i] PS3[i] ], {T, 
                  273}][[1, 2]], 
          y1[i] = x1[i] PS1[i] /P, 
y2[i] = 
            x2[i] PS2[i] /P, 
y3[i] = x3[i] PS3[i] /P, 
          x1[i + 1] = x1[i] + h (x1[i] - y1[i]), 
          x2[i + 1] = x2[i] + h (x2[i] - y2[i]), 
          x3[i + 1] = x3[i] + h (x3[i] - y3[i]), i++}];

t6 = Table[{x1[i], x2[i]}, {i, 1, 200}];

plt6 = ListPlot[t6, PlotRange -> {{0, 1}, {0, 1}}];

h = 0.1; x1[0] = 0.1; x2[0] = 0.85; x3[0] = 0.05; i = 0;

While[
        i < 200, {T = FindRoot[
                P == 
                  Evaluate[
                    x1[i] PS1[i] + x2[i] PS2[i] + x3[i] PS3[i] ], {T, 
                  273}][[1, 2]], 
          y1[i] = x1[i] PS1[i] /P, 
y2[i] = 
            x2[i] PS2[i] /P, 
y3[i] = x3[i] PS3[i] /P, 
          x1[i + 1] = x1[i] - h (x1[i] - y1[i]), 
          x2[i + 1] = x2[i] - h (x2[i] - y2[i]), 
          x3[i + 1] = x3[i] - h (x3[i] - y3[i]), i++}];

t7 = Table[{x1[i], x2[i]}, {i, 1, 200}];

plt7 = ListPlot[t7, PlotRange -> {{0, 1}, {0, 1}}];

h = 0.1; x1[0] = 0.1; x2[0] = 0.85; x3[0] = 0.05; i = 0;

While[
        i < 200, {T = FindRoot[
                P == 
                  Evaluate[
                    x1[i] PS1[i] + x2[i] PS2[i] + x3[i] PS3[i] ], {T, 
                  273}][[1, 2]], 
          y1[i] = x1[i] PS1[i] /P, 
y2[i] = 
            x2[i] PS2[i] /P, 
y3[i] = x3[i] PS3[i] /P, 
          x1[i + 1] = x1[i] + h (x1[i] - y1[i]), 
          x2[i + 1] = x2[i] + h (x2[i] - y2[i]), 
          x3[i + 1] = x3[i] + h (x3[i] - y3[i]), i++}];

t8 = Table[{x1[i], x2[i]}, {i, 1, 200}];

plt8 = ListPlot[t8, PlotRange -> {{0, 1}, {0, 1}}];

Show[plt1, plt2, plt3, plt4, plt5, plt6, plt7, plt8, 
      Graphics[{Text[C, {0, -0.1}], Text[A, {1.0, -0.1}], 
          Text[B, {0, 1.1}], RGBColor[0, 0, 1]}], AspectRatio -> 1, 
      PlotRange -> {{-0.2, 1.2}, {-0.2, 1.2}}, 
      Epilog -> {RGBColor[1, 0, 0], Line[{{1, 0}, {0, 1}}]}]

Residue Curve Map When Damkoler number (Da) equals 0.25
(i.e. slow reaction case).

Keq = 2; Da = 0.25;

h = 0.1; x1[0] = 0.2; x2[0] = 0.2; x3[0] = 0.6; i = 0;

While[
        i < 200, {T = FindRoot[
                P == 
                  Evaluate[
                    x1[i] PS1[i] + x2[i] PS2[i] + x3[i] PS3[i] ], {T, 
                  273}][[1, 2]], 
          y1[i] = x1[i] PS1[i] /P, 
y2[i] = 
            x2[i] PS2[i] /P, 
y3[i] = x3[i] PS3[i] /P, 
          x1[i + 1] = 
            x1[i] - h ((x1[i] - y1[i]) + 
                    Da (-1 + x1[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), 
          x2[i + 1] = 
            x2[i] - h ((x2[i] - y2[i]) + 
                    Da (-1 + x2[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), 
          x3[i + 1] = 
            x3[i] - h ((x3[i] - y3[i]) + 
                    Da (1 + x3[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), i++}];

t1 = Table[{x1[i], x2[i]}, {i, 1, 200}];

t1bis = Table[{x1[i], x2[i]}, {i, 1, 51}];

plt1 = ListPlot[t1, PlotRange -> {{0, 1}, {0, 1}}];

plt1bis = 
        ListPlot[t1bis, PlotRange -> {{0, 1}, {0, 1}}];

h = 0.1; x1[0] = 0.2; x2[0] = 0.2; x3[0] = 0.6; i = 0;

While[
        i < 200, {T = FindRoot[
                P == 
                  Evaluate[
                    x1[i] PS1[i] + x2[i] PS2[i] + x3[i] PS3[i] ], {T, 
                  273}][[1, 2]], 
          y1[i] = x1[i] PS1[i] /P, 
y2[i] = 
            x2[i] PS2[i] /P, 
y3[i] = x3[i] PS3[i] /P, 
          x1[i + 1] = 
            x1[i] + h ((x1[i] - y1[i]) + 
                    Da (-1 + x1[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), 
          x2[i + 1] = 
            x2[i] + h ((x2[i] - y2[i]) + 
                    Da (-1 + x2[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), 
          x3[i + 1] = 
            x3[i] + h ((x3[i] - y3[i]) + 
                    Da (1 + x3[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), i++}];

t2 = Table[{x1[i], x2[i]}, {i, 1, 200}];

plt2 = ListPlot[t2, PlotRange -> {{0, 1}, {0, 1}}];

h = 0.1; x1[0] = 0.85; x2[0] = 0.1; x3[0] = 0.05; i = 0;

While[
        i < 200, {T = FindRoot[
                P == 
                  Evaluate[
                    x1[i] PS1[i] + x2[i] PS2[i] + x3[i] PS3[i] ], {T, 
                  273}][[1, 2]], 
          y1[i] = x1[i] PS1[i] /P, 
y2[i] = 
            x2[i] PS2[i] /P, 
y3[i] = x3[i] PS3[i] /P, 
          x1[i + 1] = 
            x1[i] - h ((x1[i] - y1[i]) + 
                    Da (-1 + x1[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), 
          x2[i + 1] = 
            x2[i] - h ((x2[i] - y2[i]) + 
                    Da (-1 + x2[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), 
          x3[i + 1] = 
            x3[i] - h ((x3[i] - y3[i]) + 
                    Da (1 + x3[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), i++}];

t3 = Table[{x1[i], x2[i]}, {i, 1, 200}];

t3bis = Table[{x1[i], x2[i]}, {i, 1, 14}];

plt3 = ListPlot[t3, PlotRange -> {{0, 1}, {0, 1}}];

plt3bis = 
        ListPlot[t3bis, PlotRange -> {{0, 1}, {0, 1}}];

h = 0.1; x1[0] = 0.85; x2[0] = 0.1; x3[0] = 0.05; i = 0;

While[
        i < 200, {T = FindRoot[
                P == 
                  Evaluate[
                    x1[i] PS1[i] + x2[i] PS2[i] + x3[i] PS3[i] ], {T, 
                  273}][[1, 2]], 
          y1[i] = x1[i] PS1[i] /P, 
y2[i] = 
            x2[i] PS2[i] /P, 
y3[i] = x3[i] PS3[i] /P, 
          x1[i + 1] = 
            x1[i] + h ((x1[i] - y1[i]) + 
                    Da (-1 + x1[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), 
          x2[i + 1] = 
            x2[i] + h ((x2[i] - y2[i]) + 
                    Da (-1 + x2[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), 
          x3[i + 1] = 
            x3[i] + h ((x3[i] - y3[i]) + 
                    Da (1 + x3[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), i++}];

t4 = Table[{x1[i], x2[i]}, {i, 1, 200}];

plt4 = ListPlot[t4, PlotRange -> {{0, 1}, {0, 1}}];

h = 0.1; x1[0] = 0.1; x2[0] = 0.8; x3[0] = 0.1; i = 0;

While[
        i < 200, {T = FindRoot[
                P == 
                  Evaluate[
                    x1[i] PS1[i] + x2[i] PS2[i] + x3[i] PS3[i] ], {T, 
                  273}][[1, 2]], 
          y1[i] = x1[i] PS1[i] /P, 
y2[i] = 
            x2[i] PS2[i] /P, 
y3[i] = x3[i] PS3[i] /P, 
          x1[i + 1] = 
            x1[i] - h ((x1[i] - y1[i]) + 
                    Da (-1 + x1[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), 
          x2[i + 1] = 
            x2[i] - h ((x2[i] - y2[i]) + 
                    Da (-1 + x2[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), 
          x3[i + 1] = 
            x3[i] - h ((x3[i] - y3[i]) + 
                    Da (1 + x3[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), i++}];

t5 = Table[{x1[i], x2[i]}, {i, 1, 200}];

t5bis = Table[{x1[i], x2[i]}, {i, 1, 16}];

plt5 = ListPlot[t5, PlotRange -> {{0, 1}, {0, 1}}];

plt5bis = 
        ListPlot[t5bis, PlotRange -> {{0, 1}, {0, 1}}];

h = 0.1; x1[0] = 0.1; x2[0] = 0.8; x3[0] = 0.1; i = 0;

While[
        i < 200, {T = FindRoot[
                P == 
                  Evaluate[
                    x1[i] PS1[i] + x2[i] PS2[i] + x3[i] PS3[i] ], {T, 
                  273}][[1, 2]], 
          y1[i] = x1[i] PS1[i] /P, 
y2[i] = 
            x2[i] PS2[i] /P, 
y3[i] = x3[i] PS3[i] /P, 
          x1[i + 1] = 
            x1[i] + h ((x1[i] - y1[i]) + 
                    Da (-1 + x1[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), 
          x2[i + 1] = 
            x2[i] + h ((x2[i] - y2[i]) + 
                    Da (-1 + x2[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), 
          x3[i + 1] = 
            x3[i] + h ((x3[i] - y3[i]) + 
                    Da (1 + x3[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), i++}];

t6 = Table[{x1[i], x2[i]}, {i, 1, 200}];

plt6 = ListPlot[t6, PlotRange -> {{0, 1}, {0, 1}}];

h = 0.1; x1[0] = 0.4; x2[0] = 0.4; x3[0] = 0.2; i = 0;

While[
        i < 200, {T = FindRoot[
                P == 
                  Evaluate[
                    x1[i] PS1[i] + x2[i] PS2[i] + x3[i] PS3[i] ], {T, 
                  273}][[1, 2]], 
          y1[i] = x1[i] PS1[i] /P, 
y2[i] = 
            x2[i] PS2[i] /P, 
y3[i] = x3[i] PS3[i] /P, 
          x1[i + 1] = 
            x1[i] - h ((x1[i] - y1[i]) + 
                    Da (-1 + x1[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), 
          x2[i + 1] = 
            x2[i] - h ((x2[i] - y2[i]) + 
                    Da (-1 + x2[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), 
          x3[i + 1] = 
            x3[i] - h ((x3[i] - y3[i]) + 
                    Da (1 + x3[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), i++}];

t7 = Table[{x1[i], x2[i]}, {i, 1, 200}];

t7bis = Table[{x1[i], x2[i]}, {i, 1, 19}];

plt7 = ListPlot[t7, PlotRange -> {{0, 1}, {0, 1}}];

plt7bis = 
        ListPlot[t7bis, PlotRange -> {{0, 1}, {0, 1}}];

h = 0.1; x1[0] = 0.4; x2[0] = 0.4; x3[0] = 0.2; i = 0;

While[
        i < 200, {T = FindRoot[
                P == 
                  Evaluate[
                    x1[i] PS1[i] + x2[i] PS2[i] + x3[i] PS3[i] ], {T, 
                  273}][[1, 2]], 
          y1[i] = x1[i] PS1[i] /P, 
y2[i] = 
            x2[i] PS2[i] /P, 
y3[i] = x3[i] PS3[i] /P, 
          x1[i + 1] = 
            x1[i] + h ((x1[i] - y1[i]) + 
                    Da (-1 + x1[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), 
          x2[i + 1] = 
            x2[i] + h ((x2[i] - y2[i]) + 
                    Da (-1 + x2[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), 
          x3[i + 1] = 
            x3[i] + h ((x3[i] - y3[i]) + 
                    Da (1 + x3[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), i++}];

t8 = Table[{x1[i], x2[i]}, {i, 1, 200}];

plt8 = ListPlot[t8, PlotRange -> {{0, 1}, {0, 1}}];

h = 0.1; x1[0] = 0.0; x2[0] = 0.0; x3[0] = 1; i = 0;

FindRoot[
        P == 
          Evaluate[x1[i] PS1[i] + x2[i] PS2[i] + x3[i] PS3[i] ], {T, 
          350}][[1, 2]]

While[
        i < 200, {T = FindRoot[
                P == 
                  Evaluate[
                    x1[i] PS1[i] + x2[i] PS2[i] + x3[i] PS3[i] ], {T, 
                  273}][[1, 2]], 
          y1[i] = x1[i] PS1[i] /P, 
y2[i] = 
            x2[i] PS2[i] /P, 
y3[i] = x3[i] PS3[i] /P, 
          x1[i + 1] = 
            x1[i] + h ((x1[i] - y1[i]) + 
                    Da (-1 + x1[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), 
          x2[i + 1] = 
            x2[i] + h ((x2[i] - y2[i]) + 
                    Da (-1 + x2[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), 
          x3[i + 1] = 
            x3[i] + h ((x3[i] - y3[i]) + 
                    Da (1 + x3[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), i++}];

t9 = Table[{x1[i], x2[i]}, {i, 1, 200}];

plt9 = ListPlot[t9, PlotRange -> {{0, 1}, {0, 1}}];

h = 0.1; x1[0] = 0.0; x2[0] = 0.85; x3[0] = 0.15; i = 0;

FindRoot[
        P == 
          Evaluate[x1[i] PS1[i] + x2[i] PS2[i] + x3[i] PS3[i] ], {T, 
          350}][[1, 2]]

While[
        i < 200, {T = FindRoot[
                P == 
                  Evaluate[
                    x1[i] PS1[i] + x2[i] PS2[i] + x3[i] PS3[i] ], {T, 
                  273}][[1, 2]], 
          y1[i] = x1[i] PS1[i] /P, 
y2[i] = 
            x2[i] PS2[i] /P, 
y3[i] = x3[i] PS3[i] /P, 
          x1[i + 1] = 
            x1[i] + h ((x1[i] - y1[i]) + 
                    Da (-1 + x1[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), 
          x2[i + 1] = 
            x2[i] + h ((x2[i] - y2[i]) + 
                    Da (-1 + x2[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), 
          x3[i + 1] = 
            x3[i] + h ((x3[i] - y3[i]) + 
                    Da (1 + x3[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), i++}];

t10 = Table[{x1[i], x2[i]}, {i, 1, 200}];

plt10 = 
        ListPlot[t10, PlotRange -> {{0, 1}, {0, 1}}];

Show[plt1, plt2, plt4, plt3, plt5, plt6, plt8, plt7, plt9, plt10, 
      Graphics[{Text[C, {0, -0.1}], Text[A, {1.0, -0.1}], 
          Text[B, {0, 1.1}], RGBColor[0, 0, 1]}], AspectRatio -> 1, 
      PlotRange -> {{-0.2, 1.2}, {-0.2, 1.2}}, 
      Epilog -> {RGBColor[1, 0, 0], Line[{{1, 0}, {0, 1}}]}]

Show[plt1bis, plt2, plt4, plt3bis, plt5bis, plt6, plt8, plt7bis, plt9, 
      plt10, Graphics[{Text[C, {0, -0.1}], Text[A, {1.0, -0.1}], 
          Text[B, {0, 1.1}], RGBColor[0, 0, 1]}], AspectRatio -> 1, 
      PlotRange -> {{-0.2, 1.2}, {-0.2, 1.2}}, 
      Epilog -> {RGBColor[1, 0, 0], Line[{{1, 0}, {0, 1}}]}]

Residue Curve Map When Damkoler number (Da) equals 1.0
(i.e. intermediate reaction case).

Keq = 2; Da = 1;

h = 0.1; x1[0] = 0.2; x2[0] = 0.2; x3[0] = 0.6; i = 0;

While[
        i < 9, {T = FindRoot[
                P == 
                  Evaluate[
                    x1[i] PS1[i] + x2[i] PS2[i] + x3[i] PS3[i] ], {T, 
                  320}][[1, 2]], 
          y1[i] = x1[i] PS1[i] /P, 
y2[i] = 
            x2[i] PS2[i] /P, 
y3[i] = x3[i] PS3[i] /P, 
          x1[i + 1] = 
            x1[i] - h ((x1[i] - y1[i]) + 
                    Da (-1 + x1[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), 
          x2[i + 1] = 
            x2[i] - h ((x2[i] - y2[i]) + 
                    Da (-1 + x2[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), 
          x3[i + 1] = 
            x3[i] - h ((x3[i] - y3[i]) + 
                    Da (1 + x3[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), i++}];

t1 = Table[{x1[i], x2[i]}, {i, 1, 9}];

t1bis = Table[{x1[i], x2[i]}, {i, 1, 51}];

plt1 = ListPlot[t1, PlotRange -> {{0, 1}, {0, 1}}];

h = 0.1; x1[0] = 0.2; x2[0] = 0.2; x3[0] = 0.6; i = 0;

While[
        i < 200, {T = FindRoot[
                P == 
                  Evaluate[
                    x1[i] PS1[i] + x2[i] PS2[i] + x3[i] PS3[i] ], {T, 
                  273}][[1, 2]], 
          y1[i] = x1[i] PS1[i] /P, 
y2[i] = 
            x2[i] PS2[i] /P, 
y3[i] = x3[i] PS3[i] /P, 
          x1[i + 1] = 
            x1[i] + h ((x1[i] - y1[i]) + 
                    Da (-1 + x1[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), 
          x2[i + 1] = 
            x2[i] + h ((x2[i] - y2[i]) + 
                    Da (-1 + x2[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), 
          x3[i + 1] = 
            x3[i] + h ((x3[i] - y3[i]) + 
                    Da (1 + x3[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), i++}];

t2 = Table[{x1[i], x2[i]}, {i, 1, 200}];

plt2 = ListPlot[t2, PlotRange -> {{0, 1}, {0, 1}}];

h = 0.1; x1[0] = 0.85; x2[0] = 0.1; x3[0] = 0.05; i = 0;

While[
        i < 200, {T = FindRoot[
                P == 
                  Evaluate[
                    x1[i] PS1[i] + x2[i] PS2[i] + x3[i] PS3[i] ], {T, 
                  273}][[1, 2]], 
          y1[i] = x1[i] PS1[i] /P, 
y2[i] = 
            x2[i] PS2[i] /P, 
y3[i] = x3[i] PS3[i] /P, 
          x1[i + 1] = 
            x1[i] - h ((x1[i] - y1[i]) + 
                    Da (-1 + x1[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), 
          x2[i + 1] = 
            x2[i] - h ((x2[i] - y2[i]) + 
                    Da (-1 + x2[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), 
          x3[i + 1] = 
            x3[i] - h ((x3[i] - y3[i]) + 
                    Da (1 + x3[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), i++}];

t3 = Table[{x1[i], x2[i]}, {i, 1, 200}];

t3bis = Table[{x1[i], x2[i]}, {i, 1, 4}];

plt3 = ListPlot[t3, PlotRange -> {{0, 1}, {0, 1}}];

plt3bis = 
        ListPlot[t3bis, PlotRange -> {{0, 1}, {0, 1}}];

h = 0.1; x1[0] = 0.85; x2[0] = 0.1; x3[0] = 0.05; i = 0;

While[
        i < 200, {T = FindRoot[
                P == 
                  Evaluate[
                    x1[i] PS1[i] + x2[i] PS2[i] + x3[i] PS3[i] ], {T, 
                  273}][[1, 2]], 
          y1[i] = x1[i] PS1[i] /P, 
y2[i] = 
            x2[i] PS2[i] /P, 
y3[i] = x3[i] PS3[i] /P, 
          x1[i + 1] = 
            x1[i] + h ((x1[i] - y1[i]) + 
                    Da (-1 + x1[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), 
          x2[i + 1] = 
            x2[i] + h ((x2[i] - y2[i]) + 
                    Da (-1 + x2[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), 
          x3[i + 1] = 
            x3[i] + h ((x3[i] - y3[i]) + 
                    Da (1 + x3[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), i++}];

t4 = Table[{x1[i], x2[i]}, {i, 1, 200}];

plt4 = ListPlot[t4, PlotRange -> {{0, 1}, {0, 1}}];

h = 0.1; x1[0] = 0.1; x2[0] = 0.8; x3[0] = 0.1; i = 0;

While[
        i < 200, {T = FindRoot[
                P == 
                  Evaluate[
                    x1[i] PS1[i] + x2[i] PS2[i] + x3[i] PS3[i] ], {T, 
                  273}][[1, 2]], 
          y1[i] = x1[i] PS1[i] /P, 
y2[i] = 
            x2[i] PS2[i] /P, 
y3[i] = x3[i] PS3[i] /P, 
          x1[i + 1] = 
            x1[i] - h ((x1[i] - y1[i]) + 
                    Da (-1 + x1[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), 
          x2[i + 1] = 
            x2[i] - h ((x2[i] - y2[i]) + 
                    Da (-1 + x2[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), 
          x3[i + 1] = 
            x3[i] - h ((x3[i] - y3[i]) + 
                    Da (1 + x3[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), i++}];

t5 = Table[{x1[i], x2[i]}, {i, 1, 200}];

t5bis = Table[{x1[i], x2[i]}, {i, 1, 7}];

plt5 = ListPlot[t5, PlotRange -> {{0, 1}, {0, 1}}];

plt5bis = 
        ListPlot[t5bis, PlotRange -> {{0, 1}, {0, 1}}];

h = 0.1; x1[0] = 0.1; x2[0] = 0.8; x3[0] = 0.1; i = 0;

While[
        i < 200, {T = FindRoot[
                P == 
                  Evaluate[
                    x1[i] PS1[i] + x2[i] PS2[i] + x3[i] PS3[i] ], {T, 
                  273}][[1, 2]], 
          y1[i] = x1[i] PS1[i] /P, 
y2[i] = 
            x2[i] PS2[i] /P, 
y3[i] = x3[i] PS3[i] /P, 
          x1[i + 1] = 
            x1[i] + h ((x1[i] - y1[i]) + 
                    Da (-1 + x1[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), 
          x2[i + 1] = 
            x2[i] + h ((x2[i] - y2[i]) + 
                    Da (-1 + x2[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), 
          x3[i + 1] = 
            x3[i] + h ((x3[i] - y3[i]) + 
                    Da (1 + x3[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), i++}];

t6 = Table[{x1[i], x2[i]}, {i, 1, 200}];

plt6 = ListPlot[t6, PlotRange -> {{0, 1}, {0, 1}}];

h = 0.1; x1[0] = 0.4; x2[0] = 0.4; x3[0] = 0.2; i = 0;

While[
        i < 200, {T = FindRoot[
                P == 
                  Evaluate[
                    x1[i] PS1[i] + x2[i] PS2[i] + x3[i] PS3[i] ], {T, 
                  273}][[1, 2]], 
          y1[i] = x1[i] PS1[i] /P, 
y2[i] = 
            x2[i] PS2[i] /P, 
y3[i] = x3[i] PS3[i] /P, 
          x1[i + 1] = 
            x1[i] - h ((x1[i] - y1[i]) + 
                    Da (-1 + x1[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), 
          x2[i + 1] = 
            x2[i] - h ((x2[i] - y2[i]) + 
                    Da (-1 + x2[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), 
          x3[i + 1] = 
            x3[i] - h ((x3[i] - y3[i]) + 
                    Da (1 + x3[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), i++}];

t7 = Table[{x1[i], x2[i]}, {i, 1, 200}];

t7bis = Table[{x1[i], x2[i]}, {i, 1, 8}];

plt7 = ListPlot[t7, PlotRange -> {{0, 1}, {0, 1}}];

plt7bis = 
        ListPlot[t7bis, PlotRange -> {{0, 1}, {0, 1}}];

h = 0.1; x1[0] = 0.4; x2[0] = 0.4; x3[0] = 0.2; i = 0;

While[
        i < 200, {T = FindRoot[
                P == 
                  Evaluate[
                    x1[i] PS1[i] + x2[i] PS2[i] + x3[i] PS3[i] ], {T, 
                  273}][[1, 2]], 
          y1[i] = x1[i] PS1[i] /P, 
y2[i] = 
            x2[i] PS2[i] /P, 
y3[i] = x3[i] PS3[i] /P, 
          x1[i + 1] = 
            x1[i] + h ((x1[i] - y1[i]) + 
                    Da (-1 + x1[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), 
          x2[i + 1] = 
            x2[i] + h ((x2[i] - y2[i]) + 
                    Da (-1 + x2[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), 
          x3[i + 1] = 
            x3[i] + h ((x3[i] - y3[i]) + 
                    Da (1 + x3[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), i++}];

t8 = Table[{x1[i], x2[i]}, {i, 1, 200}];

plt8 = ListPlot[t8, PlotRange -> {{0, 1}, {0, 1}}];

h = 0.1; x1[0] = 0.0; x2[0] = 0.0; x3[0] = 1; i = 0;

While[
        i < 200, {T = FindRoot[
                P == 
                  Evaluate[
                    x1[i] PS1[i] + x2[i] PS2[i] + x3[i] PS3[i] ], {T, 
                  273}][[1, 2]], 
          y1[i] = x1[i] PS1[i] /P, 
y2[i] = 
            x2[i] PS2[i] /P, 
y3[i] = x3[i] PS3[i] /P, 
          x1[i + 1] = 
            x1[i] + h ((x1[i] - y1[i]) + 
                    Da (-1 + x1[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), 
          x2[i + 1] = 
            x2[i] + h ((x2[i] - y2[i]) + 
                    Da (-1 + x2[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), 
          x3[i + 1] = 
            x3[i] + h ((x3[i] - y3[i]) + 
                    Da (1 + x3[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), i++}];

t9 = Table[{x1[i], x2[i]}, {i, 1, 200}];

plt9 = ListPlot[t9, PlotRange -> {{0, 1}, {0, 1}}];

h = 0.1; x1[0] = 0.0; x2[0] = 0.85; x3[0] = 0.15; i = 0;

While[
        i < 200, {T = FindRoot[
                P == 
                  Evaluate[
                    x1[i] PS1[i] + x2[i] PS2[i] + x3[i] PS3[i] ], {T, 
                  273}][[1, 2]], 
          y1[i] = x1[i] PS1[i] /P, 
y2[i] = 
            x2[i] PS2[i] /P, 
y3[i] = x3[i] PS3[i] /P, 
          x1[i + 1] = 
            x1[i] + h ((x1[i] - y1[i]) + 
                    Da (-1 + x1[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), 
          x2[i + 1] = 
            x2[i] + h ((x2[i] - y2[i]) + 
                    Da (-1 + x2[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), 
          x3[i + 1] = 
            x3[i] + h ((x3[i] - y3[i]) + 
                    Da (1 + x3[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), i++}];

t10 = Table[{x1[i], x2[i]}, {i, 1, 200}];

plt10 = 
        ListPlot[t10, PlotRange -> {{0, 1}, {0, 1}}];

Show[plt1, plt2, plt4, plt3bis, plt5bis, plt6, plt8, plt7bis, plt9, 
      plt10, Graphics[{Text[C, {0, -0.1}], Text[A, {1.0, -0.1}], 
          Text[B, {0, 1.1}], RGBColor[0, 0, 1]}], AspectRatio -> 1, 
      PlotRange -> {{-0.2, 1.2}, {-0.2, 1.2}}, 
      Epilog -> {RGBColor[1, 0, 0], Line[{{1, 0}, {0, 1}}]}]

Reaction Equilibrium Curve for an ideal reacting liquid mixture (A+B=C) at Keq=2.0 .

Keq = 2; Da = 100; eps = 0.01;

i = 0; Keq = 2; h = 0.1;

While[
      i < 11, {x1[i] = 
          Evaluate[0 + h i], 
x2[
            i] = FindRoot[{x1[i] x2 - x3/Keq == 0, 
                x1[i] + x2 + x3 == 1}, {x2, 0.2}, {x3, 0.2}][[1, 
              2]], x3[
            i] = FindRoot[{x1[i] x2 - x3/Keq == 0, 
                x1[i] + x2 + x3 == 1}, {x2, 0.2}, {x3, 0.2}][[1, 
              2]], i++}]

tble = Table[{x1[i], x2[i]}, {i, 0, 10}];

Eq = 
      ListPlot[tble, PlotRange -> {{0, 1}, {0, 1}}, 
        AspectRatio -> 1, PlotJoined -> True]

Residue Curve Map When Damkoler number (Da) equals 100
(i.e. fast reaction case). Green dot in map is the singular point.

h = 0.01; x1[0] = 0.2; x2[0] = 0.2; x3[0] = 0.6; i = 0;

While[
        i < 80, {T = FindRoot[
                P == 
                  Evaluate[
                    x1[i] PS1[i] + x2[i] PS2[i] + x3[i] PS3[i] ], {T, 
                  320}][[1, 2]], 
          y1[i] = x1[i] PS1[i] /P, 
y2[i] = 
            x2[i] PS2[i] /P, 
y3[i] = x3[i] PS3[i] /P, 
          x1[i + 1] = 
            x1[i] - h (eps (x1[i] - y1[i]) +  (-1 + 
                          x1[i]) (x1[i] x2[i] - x3[i]/Keq)), 
          x2[i + 1] = 
            x2[i] - h (eps (x2[i] - y2[i]) +  (-1 + 
                          x2[i]) (x1[i] x2[i] - x3[i]/Keq)), 
          x3[i + 1] = 
            x3[i] - h (eps (x3[i] - y3[i]) + (1 + 
                          x3[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), i++}];

t1 = Table[{x1[i], x2[i]}, {i, 1, 61}];

plt1 = ListPlot[t1, PlotRange -> {{0, 1}, {0, 1}}];

h = 0.1; x1[0] = 0.2; x2[0] = 0.2; x3[0] = 0.6; i = 0;

While[
        i < 30000, {T = FindRoot[
                P == 
                  Evaluate[
                    x1[i] PS1[i] + x2[i] PS2[i] + x3[i] PS3[i] ], {T, 
                  273}][[1, 2]], 
          y1[i] = x1[i] PS1[i] /P, 
y2[i] = 
            x2[i] PS2[i] /P, 
y3[i] = x3[i] PS3[i] /P, 
          x1[i + 1] = 
            x1[i] + h (eps (x1[i] - y1[i]) + (-1 + 
                          x1[i]) (x1[i] x2[i] - x3[i]/Keq)), 
          x2[i + 1] = 
            x2[i] + h (eps (x2[i] - y2[i]) + (-1 + 
                          x2[i]) (x1[i] x2[i] - x3[i]/Keq)), 
          x3[i + 1] = 
            x3[i] + h (eps (x3[i] - y3[i]) + (1 + 
                          x3[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), i++}];

t2 = Table[{x1[i], x2[i]}, {i, 1, 30000}];

plt2 = ListPlot[t2, PlotRange -> {{0, 1}, {0, 1}}];

h = 0.01; x1[0] = 0.85; x2[0] = 0.1; x3[0] = 0.05; i = 0;

While[
        i < 200, {T = FindRoot[
                P == 
                  Evaluate[
                    x1[i] PS1[i] + x2[i] PS2[i] + x3[i] PS3[i] ], {T, 
                  273}][[1, 2]], 
          y1[i] = x1[i] PS1[i] /P, 
y2[i] = 
            x2[i] PS2[i] /P, 
y3[i] = x3[i] PS3[i] /P, 
          x1[i + 1] = 
            x1[i] - h (eps (x1[i] - y1[i]) +  (-1 + 
                          x1[i]) (x1[i] x2[i] - x3[i]/Keq)), 
          x2[i + 1] = 
            x2[i] - h (eps (x2[i] - y2[i]) +  (-1 + 
                          x2[i]) (x1[i] x2[i] - x3[i]/Keq)), 
          x3[i + 1] = 
            x3[i] - h (eps (x3[i] - y3[i]) + (1 + 
                          x3[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), i++}];

t3 = Table[{x1[i], x2[i]}, {i, 1, 200}];

t3bis = Table[{x1[i], x2[i]}, {i, 1, 55}];

plt3 = ListPlot[t3, PlotRange -> {{0, 1}, {0, 1}}];

plt3bis = 
        ListPlot[t3bis, PlotRange -> {{0, 1}, {0, 1}}];

h = 0.1; x1[0] = 0.85; x2[0] = 0.1; x3[0] = 0.05; i = 0;

While[
        i < 20000, {T = FindRoot[
                P == 
                  Evaluate[
                    x1[i] PS1[i] + x2[i] PS2[i] + x3[i] PS3[i] ], {T, 
                  273}][[1, 2]], 
          y1[i] = x1[i] PS1[i] /P, 
y2[i] = 
            x2[i] PS2[i] /P, 
y3[i] = x3[i] PS3[i] /P, 
          x1[i + 1] = 
            x1[i] + h (eps (x1[i] - y1[i]) +  (-1 + 
                          x1[i]) (x1[i] x2[i] - x3[i]/Keq)), 
          x2[i + 1] = 
            x2[i] + h (eps (x2[i] - y2[i]) +  (-1 + 
                          x2[i]) (x1[i] x2[i] - x3[i]/Keq)), 
          x3[i + 1] = 
            x3[i] + h (eps (x3[i] - y3[i]) +  (1 + 
                          x3[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), i++}];

t4 = Table[{x1[i], x2[i]}, {i, 1, 20000}];

plt4 = ListPlot[t4, PlotRange -> {{0, 1}, {0, 1}}];

h = 0.01; x1[0] = 0.1; x2[0] = 0.8; x3[0] = 0.1; i = 0;

While[
        i < 200, {T = FindRoot[
                P == 
                  Evaluate[
                    x1[i] PS1[i] + x2[i] PS2[i] + x3[i] PS3[i] ], {T, 
                  273}][[1, 2]], 
          y1[i] = x1[i] PS1[i] /P, 
y2[i] = 
            x2[i] PS2[i] /P, 
y3[i] = x3[i] PS3[i] /P, 
          x1[i + 1] = 
            x1[i] - h (eps (x1[i] - y1[i]) + (-1 + 
                          x1[i]) (x1[i] x2[i] - x3[i]/Keq)), 
          x2[i + 1] = 
            x2[i] - h (eps (x2[i] - y2[i]) +  (-1 + 
                          x2[i]) (x1[i] x2[i] - x3[i]/Keq)), 
          x3[i + 1] = 
            x3[i] - h (eps (x3[i] - y3[i]) +  (1 + 
                          x3[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), i++}];

t5 = Table[{x1[i], x2[i]}, {i, 1, 200}];

t5bis = Table[{x1[i], x2[i]}, {i, 1, 126}];

plt5 = ListPlot[t5, PlotRange -> {{0, 1}, {0, 1}}];

plt5bis = 
        ListPlot[t5bis, PlotRange -> {{0, 1}, {0, 1}}];

h = 0.1; x1[0] = 0.1; x2[0] = 0.8; x3[0] = 0.1; i = 0;

While[
        i < 10000, {T = FindRoot[
                P == 
                  Evaluate[
                    x1[i] PS1[i] + x2[i] PS2[i] + x3[i] PS3[i] ], {T, 
                  273}][[1, 2]], 
          y1[i] = x1[i] PS1[i] /P, 
y2[i] = 
            x2[i] PS2[i] /P, 
y3[i] = x3[i] PS3[i] /P, 
          x1[i + 1] = 
            x1[i] + h (eps (x1[i] - y1[i]) +  (-1 + 
                          x1[i]) (x1[i] x2[i] - x3[i]/Keq)), 
          x2[i + 1] = 
            x2[i] + h (eps (x2[i] - y2[i]) +  (-1 + 
                          x2[i]) (x1[i] x2[i] - x3[i]/Keq)), 
          x3[i + 1] = 
            x3[i] + h (eps (x3[i] - y3[i]) +  (1 + 
                          x3[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), i++}];

t6 = Table[{x1[i], x2[i]}, {i, 1, 10000}];

plt6 = ListPlot[t6, PlotRange -> {{0, 1}, {0, 1}}];

h = 0.1; x1[0] = 0.4; x2[0] = 0.4; x3[0] = 0.2; i = 0;

While[
        i < 200, {T = FindRoot[
                P == 
                  Evaluate[
                    x1[i] PS1[i] + x2[i] PS2[i] + x3[i] PS3[i] ], {T, 
                  273}][[1, 2]], 
          y1[i] = x1[i] PS1[i] /P, 
y2[i] = 
            x2[i] PS2[i] /P, 
y3[i] = x3[i] PS3[i] /P, 
          x1[i + 1] = 
            x1[i] - h (eps (x1[i] - y1[i]) + (-1 + 
                          x1[i]) (x1[i] x2[i] - x3[i]/Keq)), 
          x2[i + 1] = 
            x2[i] - h (eps (x2[i] - y2[i]) +  (-1 + 
                          x2[i]) (x1[i] x2[i] - x3[i]/Keq)), 
          x3[i + 1] = 
            x3[i] - h (eps (x3[i] - y3[i]) +  (1 + 
                          x3[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), i++}];

t7 = Table[{x1[i], x2[i]}, {i, 1, 200}];

t7bis = Table[{x1[i], x2[i]}, {i, 1, 14}];

plt7 = ListPlot[t7, PlotRange -> {{0, 1}, {0, 1}}];

plt7bis = 
        ListPlot[t7bis, PlotRange -> {{0, 1}, {0, 1}}];

h = 0.1; x1[0] = 0.4; x2[0] = 0.4; x3[0] = 0.2; i = 0;

While[
        i < 20000, {T = FindRoot[
                P == 
                  Evaluate[
                    x1[i] PS1[i] + x2[i] PS2[i] + x3[i] PS3[i] ], {T, 
                  273}][[1, 2]], 
          y1[i] = x1[i] PS1[i] /P, 
y2[i] = 
            x2[i] PS2[i] /P, 
y3[i] = x3[i] PS3[i] /P, 
          x1[i + 1] = 
            x1[i] + h (eps (x1[i] - y1[i]) +  (-1 + 
                          x1[i]) (x1[i] x2[i] - x3[i]/Keq)), 
          x2[i + 1] = 
            x2[i] + h (eps (x2[i] - y2[i]) +  (-1 + 
                          x2[i]) (x1[i] x2[i] - x3[i]/Keq)), 
          x3[i + 1] = 
            x3[i] + h (eps (x3[i] - y3[i]) +  (1 + 
                          x3[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), i++}];

t8 = Table[{x1[i], x2[i]}, {i, 1, 20000}];

plt8 = ListPlot[t8, PlotRange -> {{0, 1}, {0, 1}}];

h = 0.1; x1[0] = 0.0; x2[0] = 0.0; x3[0] = 1; i = 0;

While[
        i < 30000, {T = FindRoot[
                P == 
                  Evaluate[
                    x1[i] PS1[i] + x2[i] PS2[i] + x3[i] PS3[i] ], {T, 
                  273}][[1, 2]], 
          y1[i] = x1[i] PS1[i] /P, 
y2[i] = 
            x2[i] PS2[i] /P, 
y3[i] = x3[i] PS3[i] /P, 
          x1[i + 1] = 
            x1[i] + h (eps (x1[i] - y1[i]) +  (-1 + 
                          x1[i]) (x1[i] x2[i] - x3[i]/Keq)), 
          x2[i + 1] = 
            x2[i] + h (eps (x2[i] - y2[i]) +  (-1 + 
                          x2[i]) (x1[i] x2[i] - x3[i]/Keq)), 
          x3[i + 1] = 
            x3[i] + h (eps (x3[i] - y3[i]) +  (1 + 
                          x3[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), i++}];

t9 = Table[{x1[i], x2[i]}, {i, 1, 30000}];

plt9 = ListPlot[t9, PlotRange -> {{0, 1}, {0, 1}}];

h = 0.1; x1[0] = 0.0; x2[0] = 0.85; x3[0] = 0.15; i = 0;

While[
        i < 20000, {T = FindRoot[
                P == 
                  Evaluate[
                    x1[i] PS1[i] + x2[i] PS2[i] + x3[i] PS3[i] ], {T, 
                  320}][[1, 2]], 
          y1[i] = x1[i] PS1[i] /P, 
y2[i] = 
            x2[i] PS2[i] /P, 
y3[i] = x3[i] PS3[i] /P, 
          x1[i + 1] = 
            x1[i] + h (eps (x1[i] - y1[i]) +  (-1 + 
                          x1[i]) (x1[i] x2[i] - x3[i]/Keq)), 
          x2[i + 1] = 
            x2[i] + h (eps (x2[i] - y2[i]) +  (-1 + 
                          x2[i]) (x1[i] x2[i] - x3[i]/Keq)), 
          x3[i + 1] = 
            x3[i] + h (eps (x3[i] - y3[i]) +  (1 + 
                          x3[i]) (x1[i] x2[i] - 
                          x3[i]/Keq)), i++}];

t10 = Table[{x1[i], x2[i]}, {i, 1, 20000}];

plt10 = 
        ListPlot[t10, PlotRange -> {{0, 1}, {0, 1}}];

Show[plt3bis, plt1, plt2, plt4, plt6, plt5bis, plt7bis, plt8, plt9, 
        plt10, Eq, 
        Graphics[{PointSize[0.04], RGBColor[0, 1, 0], 
            Point[{0.11396996171408827, 0.7200396519803999}]}], 
        AspectRatio -> 1, 
        Epilog -> {RGBColor[1, 0, 0], 
            Line[{{1, 0}, {0, 1}}]}];

Determination of Singular Points for Da=0.25, 1, 100 .

T =.

A1 = 23.1962; B1 = -3502.77; C1 = 0.0;
A2 = 22.6867; B2 = -3502.77; C2 = 0.0;
A3 = 21.5872; B3 = -3502.77; C3 = 0.0;

PS1 = Exp[A1 + B1/(C1 + T)];
PS2 = Exp[A2 + B2/(C2 + T)];
PS3 = Exp[A3 + B3/(C3 + T)];

Da = 100;

FindRoot[{P == 
          x1 PS1 + 
            x2 PS2 + (1 - x1 - x2) PS3, (P x1 - 
                x1 PS1) (-1 + x2) ==  (P x2 - 
                x2 PS2) (-1 + x1), 
(P x2 - 
              x2 PS2) == -  Da P (x1 x2 - (1 - x1 - 
                      x2)/Keq) (-1 + x2)}, {x1, 0.09, 0.12}, {x2, 
        0.7, 0.75}, {T, 320, 360}, MaxIterations -> 50, 
      AccuracyGoal -> 5]

Da = 1;

FindRoot[{P == 
          x1 PS1 + 
            x2 PS2 + (1 - x1 - x2) PS3, (P x1 - 
                x1 PS1) (-1 + x2) ==  (P x2 - 
                x2 PS2) (-1 + x1), 
(P x2 - 
              x2 PS2) == -  Da P (x1 x2 - (1 - x1 - 
                      x2)/Keq) (-1 + x2)}, {x1, 0.1, 0.2}, {x2, 
        0.3, 0.4}, {T, 320, 360}, MaxIterations -> 50, 
      AccuracyGoal -> 5]

Da = 0.25;

FindRoot[{P == 
          x1 PS1 + 
            x2 PS2 + (1 - x1 - x2) PS3, (P x1 - 
                x1 PS1) (-1 + x2) ==  (P x2 - 
                x2 PS2) (-1 + x1), 
(P x2 - 
              x2 PS2) == -  Da P (x1 x2 - (1 - x1 - 
                      x2)/Keq) (-1 + x2)}, {x1, 0.01, 0.05}, {x2, 
        0.05, 0.2}, {T, 320, 360}, MaxIterations -> 50, 
      AccuracyGoal -> 5]