Issue plotting maximum and minimum roots of polynomial

So I'm basically writing a script that needs to find the maximum and minimum roots of a three root polynomial, for a range of input values. Then, I want to plot these maximum and minimum roots on a graph to obtain a curve.

Now, in the current script I have written it is plotting all of the MAXIMUM roots fine, but for some reason it isn't plotting the MINIMUM roots correctly. The minimum roots should generally be in the range of 40-50 on the x-axis. Here is my script:

 clc;
clear;
Tc = 126.2; 
Pc = 34; 
Omega = 0.0377; 
m = 0.37464 + 1.54226*Omega - 0.26992*Omega^2;
A = 3.7362;
B = 264.651;
C = -6.788; 
j = 1; 
R = 83.14; 
b = 0.07780*(R*Tc/Pc);
for i = 1:1:63;
T(i) = 63.2 + i; 
Tre = T(i)/Tc; 
f = A - (B/(T(i)+C)); 
P = 10.^f; 
a = 0.45724*(R*Tc)^2/Pc; 
alph = (1+m*(1-sqrt(Tre)))^2; 
PR = [P (P.*b - R.*T) (a.*alph - 3.*b.^2*P - 2.*b.*R.*T) b.*(b.^2.*P + b.*R.*T - a.*alph)];
Vsat = roots(PR); 
Vsat = Vsat(imag(Vsat)==0); 
Vliq = min(Vsat);
Vvap = max(Vsat); 
plot(Vliq,P,'.');
hold on;
plot(Vvap,P,'.');
axis([0 500 0 80]); 
  end

Attached is a picture of the curve I am getting and the curve I should be getting. All of the point on the right (max roots) are fine, whereas the ones on the left (min roots) are wrong for some reason. Any help would be great thanks.

[Adam]

Is the data correct? The results of a calculation and the results of a plot are two separate things and the first thing you should do when your plot does not look as expected is to ascertain whether it is the plotting function that is doing something unexpected or whether the data you are giving it is incorrect. You can check this easily using breakpoints and the variable editor or command line.

[Peter]

Yeah I've written another script where I'm putting in individual values and it is returning values that are much more reasonable and where they should be. I copied over the polynomial data to this script so it is identical. I tried using breakpoints and for some reason I can see it is returning negative values for the minimum, but values which are fine the for the maximum. Very confused.

Answer 1

[David Goodmanson]

Hi Peter,

I believe that the basic problem here is the three lines of code

T(i) = 63.2 + i; 
Tre = T(i)/Tc; 
f = A - (B/(T(i)+C));

rather than

T = 63.2 + i; 
Tre = T/Tc; 
f = A - (B/(T+C));

The first way makes T into an array that increases in size as the for loop progresses, giving you more and more roots as you go along. There still appears to be something not quite right with the algebra since the result does not display the critical point at T = 126.2, but at least it's one issue out of the way.