Why is my plot not showing me the right curve

I am doing a thermo project to find the heat capacities of a compound, and in reality, the CP-CV curve should only be positive. I am pretty sure that my parameters are correct. I don't understand why is the plot not smooth. Can someone check my code, if i am doing anything incorrect? Thank you!

Ttp = 148.35;
Tc = 421.9;%K from NIST
Pc = 2177490;%Pa from Chemeo
R = 8.3141;%L/(K*mol)

%Antoine Parameters and Functions
antA = 4.38423;
antB = 1257.758;
antC = -13.231;

syms T 
antP = 10^(antA-(antB/(T+antC)));%antoine pressure formula
Pfxn = matlabFunction(antP);
antT =0.7*Tc;
P = Pfxn(antT)* 10^5;
w = -1 -log10(P/Pc)
k = 0.37464 + 1.54226*w - 0.26992*w^2

sqrtAlpha = 1 +k*(1-sqrt(T/Tc));
al = sqrtAlpha^2
a = 0.45724*(R^2)*(Tc^2)*(1/Pc)*al;
b = 0.0778*R*Tc/Pc;

syms V 
PRp = ((R*T)/(V-b)) - (a/((V*(V+b)+b*(V-b))));

%setting up the temperature range
Trange = Ttp:1:(2.5*Tc);
[m,n] = size(Trange)
Vrange = []; %setting as a matrix

for i = Trange
    Vsolve = vpasolve([T == i,PRp == Pc], [T,V])
    Vsolve = real(double(Vsolve.V))
    Vrange(end+1) = Vsolve
 
end

dpdtV = diff(PRp,T);
dpdvT = diff(PRp,V);
CpCvf = T*(-(dpdtV^2)/dpdvT);
CpCvfxn = matlabFunction(CpCvf);
CpCv = CpCvfxn(Trange,Vrange)

plot(CpCv)
xlim([200 1000])
ylim([0 1300])
xlabel("Temperature (K)")
ylabel("Cp-Cv (J/mol*K)")
title("Cp-Cv vs. T (Perfluoropentane), Pr=1.0")

[darova]

Where is the error?

[Minwei Lin]

there is no error message, just that the plot generates a negative value while there should not be one

Answer 1

[Walter Roberson]

    Vsolve = vpasolve([T == i,PRp == Pc], [T,V])    

In that statement, T is a symbolic variable and i is a numeric scalar, and Pc is a numeric constant. The combination effectively ends up defining V as a cubic. However, vpasolve() ends up selecting one of the roots -- and you have no control over which root it selects.

Your use of real() around the solution in the line

Vsolve = real(double(Vsolve.V))

tells us that you are expecting that the solution might be complex valued. If so, then which of the three roots are you expecting?

The equation has two complex-valued solutions and one real-valued solution until i becomes 803.35, at which point it starts having three real-valued solutions, two of which are small and negative. vpasolve(), in picking one of the three, happens to pick one of the negative ones.

I would suggest to you that you do not want complex-valued solutions, that instead you want the maximum of whatever real-valued solutions are available. I would thus suggest,

eqn = solve(PRp==Pc,V);
for i = Trange
    Vsolve = double(subs(eqn, T, i));
    Vreal = Vsolve(imag(Vsolve) == 0);
    if ~isempty(Vreal)
        Vsolve = max(Vreal);
    else
        Vsolve = real(Vsolve(1));  %not reached in practice
    end
    Vrange(end+1) = Vsolve;
end

[Minwei Lin]

thank you for your answer!! I solved it already and it's exactly what you suggest.