Trial and error help is not giving me solution to that (still waiting ....

% Calculate Bubble Point and Dew Point SRK-EOS
% P = RT / V - b  -  alpha*a / V(V-b)
%----------------------------------------------

% liquid phase composition

y = [0.86 0.05 0.05 0.02 0.01 0.0095 0.0005]';

% vapor phase composition

x = [0.45 0.05 0.05 0.03 0.01 0.01 0.4]';

% critical pressure

Pc =[4.5947 4.8711 4.2472 3.6397 3.3688 3.1233 1.965]';

% critical temperature

Tc =[190.7389 305.5111 370.0333 425.3444 469.8889 512.7778 644.4444]';

% Accentric Factor

w = xlsread('myfile.xlsx','Sheet1','E2:E8');

%---------------------------------------------

L = 1;

R = 8.314e-3; % MPa m3 / mole.K

T = 200; % K

Tr = T./Tc; % K

P(L) = 5.0; %MPa

m = 0.480+1.574.*w -0.176.*w.^2;

al=1+ m.*(1-sqrt(Tr)).^2;

a = 0.42748.*(((R.*Tc).^2)./Pc);

b = 0.08664.*(R.*Tc./Pc);

q=0;

while q==0

      % Z3 - Z2 + (A-B-B^2)Z-AB=0

      % for gas phase

      SA = 0;

      RR = zeros(7,1);

      for i = 1:7

          for j = 1:7

              SA = SA + y(i).*y(j).*sqrt(a(i).*a(j).*al(i).*al(j));

              RR(i) = RR(i) + y(j).*sqrt(a(i).*a(j).*al(i).*al(j));

          end

      end

      SB = sum(y.*b);

      % For Liquid Phase

      SAA = 0;

      MM = zeros(7,1);

      for i = 1:7

          for j = 1:7

              SAA = SAA + x(i).*x(j)*sqrt(a(i).*a(j)*al(i).*al(j));

              MM(i) = MM(i) +x(j).*sqrt(a(i).*a(j).*al(i).*al(j));

          end

      end

      SBB = sum(x.*b);

      %for Gas Phase

      A = SA.*P(L)./(R*T).^2;

      B = SB.*P(L)./(R*T);

      ZV=roots([1 -1 A-B-B^2 -A*B]);

      Zv=max(ZV);

      Zvc=isreal(Zv);

      if Zvc == false

          Zv=abs(Zv);

          disp('error 1')

      end

      %for liquid Phase

      AA = (SAA.*P(L))./(R*T).^2;

      BB = (SBB.*P(L)./(R*T));

      Z=roots([1 -1  AA-BB-BB^2 -AA*BB]);

      Zl=min(Z);

      Zvl=isreal(Zl);

      if Zvl == false

          Zl=abs(Zl);

          disp('error 2')

      end

      phiv =exp((b.*(Zv-1)/SB) - log(Zv-B)-(A/B).*((2.*RR/SA)-

(b/SB)).*log(1 + B/Zv));

      phil =exp((b.*(Zl-1)./SBB) - log(Zl-BB)-(AA./BB).*((2.*MM/SAA)-

(b./SBB)).*log(1 + BB./Zl));

      K = phil./phiv;

      RA1= sum(y./K)-1;

      if abs(RA1) < 1e-5 && abs((y./K) - x) < 1e-5

          break

      else

          if L==1

              L = L+1;

              P(L) = 1.01.*P(L-1);

              RA1;

          else

              if L==1

                  P(L)=P(L-1)-(RA1./(RA1-RA2))*(P(L-1)-P(L-2));

                  RA2 = RA1;

              end

              %adjustment

              x2 = (1-0.5)*(y./K) + 0.5*x;

              x = x2;

          end

      end

end

disp(P(L));

% I am trying to using trial and error to get the results of P(L) , but my solution is not converging and keep going on. I need to know what can I do to resolve the issue .

Help Please