Trjtdtk

Maxwell Construction is just integral of area and finding root to find vapor pressure.

This must be exactly the same where Gibbs energy graph cross at a point i found it p = 0.81.

The problem my areadifferential doesn't work properly to return a value p=0.81;

What mistake i have done in coding!?

t = 0.95; 

p[v_] := (8*t)/(3*v - 1) - 3/v^2;

Plot[p[v], v, 0.5, 3, PlotRange -> 0, 3, 0.6, 1, AxesLabel -> V/Vc, P/Pc]


g[v_] := (-t)*Log[3*v - 1] + 0.95/(3*v - 1) - 9/(4*v); 

ParametricPlot[p[v],g[v],v,.65,2.25,AxesLabel->P/Pc,G/NKT,
PlotLabel->"Gibbs Free Energy Vs. P/Pc"]

That's Return pressure p=P/Pc=0.81 in plot which is correct

This is the error Maxwell Construction "area differential" Part

pint[v_] := (8/3)*t*Log[3*v - 1] + 3/v;

areadifferential[p0_, v1guess_, v2guess_] := 
 (v1 = FindRoot[p[v] == p0, v, v1guess][[1,2]]; 
 v2 = FindRoot[p[v] == p0, v, v2guess][[1,2]];
pint[v2] - pint[v1] - p0*(v2 - v1))

FindRoot[areadifferential[p0, 0.7, 2] == 0, p0, 0.8, 0.82]

Good job on describing your problem.

Copying your definitions from your question

p[v_] := (8 t)/(3 v - 1) - 3/v^2
g[v_] := t/(3 v - 1) - t Log[-1 + 3 v] - 9/(4 v)
pint[v_] := 8/3 t Log[-1 + 3 v] + 3/v
t = 0.95

I found the value 0.81188 to be the value that minimized your areadifferential function (shown at the bottom).

Look at the plot of P/Pc vs V/Vc and note the three points corresponding to three different V/Vc values resulting in p[v] = 0.8118.

Module[
 
 v1 = v /. FindRoot[p[v] == 0.8118, v, 0.5],
 p1,
 v2 = v /. FindRoot[p[v] == 0.8118, v, 1],
 p2,
 v3 = v /. FindRoot[p[v] == 0.8118, v, 2],
 p3
 ,
 p1 = p[v1];
 p2 = p[v2];
 p3 = p[v3];
 Show[
 Plot[p[v], v, 0.65, 2.25, PlotRange -> 0, 3, 0.6, 1, 
 AxesLabel -> "V/Vc", "P/Pc"],
 Graphics[
 
 PointSize[0.03],
 Red,
 Point[v1, p1],
 Green,
 Point[v2, p2],
 Black,
 Point[v3, p3]
 
 ]
 ]
 ]

The problem that you are experiencing is because the value of 0.7 for v1guess is just to the right of the valley so it converges on the wrong point (i.e., the green point).

Change v1guess to 0.5 and you will be fine, it will converge on the red point. Using v2guess of 2.0 works fine causing convergence on the black point.

I slightly modified the areadifferential function to use Module rather than parenthesis. One needs to constrain the input arguments to be numeric so FindRoot doesn't complain (FindFoot first tries to work with symbolic arguments which doesn't work with your problem).

areadifferential[
p0_?NumericQ,
v1guess_?NumericQ,
v2guess_?NumericQ] := Module[
 
 v1 = FindRoot[p[v] == p0, v, v1guess][[1, 2]],
 v2 = FindRoot[p[v] == p0, v, v2guess][[1, 2]]
 ,
 pint[v2] - pint[v1] - p0*(v2 - v1)
 ]

and then

FindRoot[areadifferential[p0, 0.5, 2] == 0, p0, 0.8, 0.82]
(* p0 -> 0.811879 *)