Newton Raphson apply to Van der Waal Equation

Question

I need assistance with this problem. I have two functions, where f1 corresponds to the dimensionless van der Waals equation of state, and function 2 corresponds to the dimensionless chemical potential function. Xl corresponds to the liquid mole fraction, and Xv corresponds to the vapor mole fraction. Y is equal to T/Tc. Z is equal to P/Pv and is also the dimensionless van der Waals equation of state. I have this code, and I need the Newton-Raphson method with 1st order continuation to solve this equation for values of 0.50 < y < 1.0. As an initial step, I used x = {0.89, 1.2}. I need to obtain an isotherm at the end of the code, but my code is running infinitely and don't find a convergence. In this case I need to fix y and obtain the values of xl and xv and then apply to Z to obtain the isotherm

f1 = -(3/xl^2) + 3/xv^2 + (8 *(1/(-1 + (3 *xl)) + 1/(1 - (3* xv))) *y);
f2 = (-3 (3 - 9 xl + 4 xl^2 y) + 4 xl (-1 + 3 xl) y Log[-1 + 3 xl])/(
   xl (-1 + 3 xl) (3 + Log[16])) - (-3 (3 - 9 xv + 4 xv^2 y) + 
    4 xv (-1 + 3 xv) y Log[-1 + 3 xv])/(xv (-1 + 3 xv) (3 + Log[16]));


f = {f1, f2};
J = D[f, {{xl, xv}}];
dfdy = D[f, y];
yIter = 0.999;
iterations = 100;
tolerance = 1*^-6;
step = 0.001;

size = Length[Range[0.5, 1, step]];
Z = ConstantArray[0, {2, size}];
X = ConstantArray[0, {2, size}];
xInit = {0.89, 1.20};

For[index = 1, index <= size, index++, yIter = 1 - index*step;
  For[iter = 1, iter <= iterations, iter++, 
   JEval = N[J /. {xl -> xInit[[1]], xv -> xInit[[2]], y -> yIter}];
   JInv = Inverse[JEval];
   fEval = N[f /. {xl -> xInit[[1]], xv -> xInit[[2]], y -> yIter}];
   deltaX = JInv.-fEval;
   xInit = xInit + deltaX;
   If[Norm[fEval] < tolerance, X[[All, index]] = xInit;
    Z[[All, index]] = 8*yIter/(3*xInit[[1]] - 1) - (3/xInit[[1]]^2);
    deltaX = 
     JInv.(step*
        D[dfdy /. {xl -> xInit[[1]], xv -> xInit[[2]], y -> yIter}]);
    xInit = xInit + deltaX;
    Break[];];
   If[iter == 99, Print["Did not converge"]];];];

Export["vle.svg", 
  ListLinePlot[Transpose[{X[[1]], Z[[1]]}], PlotStyle -> Black]];
Export["vle.svg", 
 ListLinePlot[Transpose[{X[[2]], Z[[2]]}], PlotStyle -> Black]]

I need to find this curve. The point on the curve denotes the critical point. The left branch of the curve represents the molar volume of coexisting liquid and the right branch represents the molar volume of the vapor. The points defining the curve are generated by solving the equilibrium equation for various temperatures. Note that are two curves (left and right side).

Answer 1

To find curve we can use this code

f1 = -(3/xl^2) + 3/xv^2 + (8*(1/(-1 + (3*xl)) + 1/(1 - (3*xv)))*y);
f2 = (-3 (3 - 9 xl + 4 xl^2 y) + 
      4 xl (-1 + 3 xl) y Log[-1 + 3 xl])/(xl (-1 + 3 xl) (3 + 
        Log[16])) - (-3 (3 - 9 xv + 4 xv^2 y) + 
      4 xv (-1 + 3 xv) y Log[-1 + 3 xv])/(xv (-1 + 3 xv) (3 + 
        Log[16]));


f = {f1, f2}; var = {xl, xv};
J = D[f, {var}];
dfdy = D[f, y];
yIter = 0.999;
iterations = 100;
tolerance = 1*^-6;
step = 0.001;

size = Length[Range[0.5, 1, step]];
Z = ConstantArray[0, {2, size}];
X = ConstantArray[0, {2, size}];
xInit = {0.89, 1.20};


For[index = 1, index <= size, index++, yIter = 1 - index*step;
  For[iter = 1, iter <= iterations, iter++, 
   JEval = J /. {xl -> xInit[[1]], xv -> xInit[[2]], y -> yIter};
   JInv = Inverse[JEval];
   fEval = f /. {xl -> xInit[[1]], xv -> xInit[[2]], y -> yIter};
   deltaX = JInv . (-fEval);
   xInit = xInit + deltaX;
   If[Norm[fEval] < tolerance, X[[All, index]] = xInit;
    Z[[All, index]] = 8*yIter/(3*xInit[[1]] - 1) - (3/xInit[[1]]^2);
    deltaX = 
     JInv . (step*
        D[dfdy /. {xl -> xInit[[1]], xv -> xInit[[2]], y -> yIter}]);
    xInit = xInit + deltaX;
    Break[];];
   If[iter == 99, Print["Did not converge"]];];];

Visualization

    {ListLinePlot[Transpose[{X[[1]], Z[[1]]}], PlotStyle -> Black, 
      PlotRange -> All], 
     ListLinePlot[Transpose[{X[[2]], Z[[2]]}], PlotStyle -> Black, 
      PlotRange -> All], ListPlot[
 Join[Transpose[{X[[1]], Z[[1]]}], Transpose[{X[[2]], Z[[2]]}]], 
 PlotRange -> All, PlotStyle -> Black, Frame -> True, FrameLabel -> {"V/Vc", "P/Pc"}]}