Q = @(v) sym(v);
Tc = Q(400);
Pc = Q(35);
R = Q(0.0821);
Z = Q(0.34);
syms Vc a b c;
%assume([b, c], 'real')
assumeAlso(b ~= 0 & c ~= 0)
X = R*Tc*(1-c)/(Vc)-(a*b*(c^2)/(Vc^2))-(R*Tc/b)*log(1-((c*b)/Vc))-Pc;
eqn1 = R*Tc*(1-c)/(Vc)-(a*b*(c^2)/(Vc^2))-(R*Tc/b)*log(1-((c*b)/Vc))-Pc == 0
eqn2 = diff(X,Vc) == 0
eqn3 = diff(X,Vc,2) == 0
eqn4 = (Pc*Vc)/(R*Tc)-Z == 0
eqns = simplify([eqn1, eqn2, eqn3, eqn4])
Vc_sol = solve(eqns(end), Vc)
eqns2 = simplify(subs(eqns(1:end-1), Vc, Vc_sol))
partial_a = solve(eqns2(2), a)
eqns3 = simplify(subs(eqns2([1 3:end]), a, partial_a))
partial_b = solve(eqns3(2), b, 'returnconditions', true)
partial_b.b
partial_b.conditions
eqns4_1 = (subs(eqns3([1 3:end]), b, partial_b.b(1)))
eqns4_2 = (subs(eqns3([1 3:end]), b, partial_b.b(2)))
sol_c_1 = vpasolve(eqns4_1, c)
sol_c_2 = vpasolve(eqns4_2, c)
full_c = sol_c_2
full_b = subs(partial_b.b(2), c, full_c)
full_a = subs(subs(partial_a, b, partial_b.b(2)), c, full_c)
full_Vc = Vc_sol
sol = [full_a, full_b, full_c, full_Vc]
subs([eqn1, eqn2, eqn3, eqn4], [a, b, c, Vc], sol)
vpa(ans)
So the solution works to within round-off error.
