for and while loops and graphing

Tc = 670.2; %in K
Pc = 55.02; %in Bar
Omega = 0.39983;
Kappa1 = -0.03471;
Kappa0 = 0.378893 + 1.4897153*Omega - 0.17131848*Omega^2 + 0.0196554*Omega^3;
R = 0.00008314; %Gas constant
b = .0778*R*Tc/Pc; %b is constant, so can leave out the loop
P = 0.0000154565:.001:Pc; %create a range of pressures from 0 to critical pressure
T = 236.65:1:Tc; %creates a range of temps from close to 0 to critical temp
Tr = T/Tc;
ratio = 1;
K = 0;
for m = 1:length(T)%goes through every tempurature from T range
 for j = 1:length(P)%tests every pressure from pressure range at a certain T to find where fugacities are equal
   while ratio < 0.0001
     Kappa = Kappa0 + Kappa1*(1 + Tr.^0.5).*(0.7-Tr);
     alpha = (1+Kappa*(1-(sqrt(T(m)/Tc))))^2;
     a = .045724*(R^2)*(Tc^2)/Pc*alpha;
     A = a*P(j)/((R^2)*(T(m)^2));
     B = b*P(j)/(R*T(m));
     y = [(1) -(1-B) (A-3*B^2-2*B) -(A*B-B^2-B^3)]; %Peng-Robinson Cubic Root
     r = roots(y); %Solves for the roots in the P-R cubic
     Zl = (min(r)); %Lowest root is the compressibility for liquid
     Zv = (max(r)); %Highest root is compressibility for vapor
     fv = ((P(j)*exp((Zv-1)-log(Zv-B)-A/(2*sqrt(2)*B)*log((Zv+(1+sqrt(2))*B)/(Zv+(1-sqrt(2))*B))))); %fugacity of vapor
     fl = ((P(j)*exp((Zl-1)-log(Zl-B)-A/(2*sqrt(2)*B)*log((Zl+(1+sqrt(2))*B)/(Zl+(1-sqrt(2))*B))))); %fugacity of liquid
     ratio = abs(fl/fv-1); %calculates how close fugacity of liquid and vapor are
     P(j) = P(j)*fl/fv
     K = K +1;
   end
 end
end

The first loop is to go through the temperature range, and the second one is for the pressure. The while loop is the criterion of stability. This is what I want the program to do with this example: for the first Temperature which is 236.65 and first pressure which is 0.0000154565. I will use those two to find the roots of a cubic equation, then calculate the fugacity of the max root and the min root. If the fugacities of the two do not equal, i will multiply the pressure 0.0000154565 with fl/fv to get a new pressure and this new pressure along with the unchanged first temperature 236.65 will be used to recalculate fugacities again until the condiditon is satisfied. Once the computations for each Temperature and Pressure are completed, I want to plot a T(xaxis) vs P(yaxis) graph. How to plot the graph? where to put the plot function?

Thanks to anyone who took the time to read my question.