MATLAB code never stops running (infinite loop)

My Matlab code stuck in infinite loop, error comes from line 157 to 242 (fourth while loop). I really don't know what is wrong. I would appreciate any help on this!

clear all
clear variables
global nos noc nof CompMW CompSigma FeedTemp FeedPres Ftotal ...
        length ir functiontemp
errormax=1e-6;
maxiterations=500;
nos=100;
noc=2;
z(1)=0.21;
z(2)=0.79;
Q(2)=1.04e-6;
selectivity=5.85;
Q(1)=Q(2)*selectivity;
FeedPres=150*76/14.7;  %psia to cmHg
PermPres=14.7*76/14.7;
FeedTemp=25+273.15;
CompMW(1)=32;  %O2 (g/mol)
CompMW(2)=28.02;  %N2
CompMW(3)=44.01;   %CO2
CompSigma(1)=3.433; %Angstroms
CompSigma(2)=3.667;
CompSigma(3)=3.996;
Ptotal=zeros(1,nos);  %Pre-allocations for faster code execution
Rtotal=zeros(1,nos);
R=zeros(nos,noc);
P=zeros(nos,noc);
Sweep=zeros(1,noc);
R_temp=zeros(nos,noc);
P_temp=zeros(nos,noc);
nof=1000; %number of fibers
length=50; %cm
or=0.04/2; %cm
now=6;
fmin=nof*2*pi*or*length*Q(1)*FeedPres;
for n=1:now
    Ftotal(n,1)=fmin;
end
PresDropSet=0.04;  %psia
PresDropMax=0.2;

x(1)=-2.3506;
x(2)=-1.33584;
x(3)=-0.43607;
x(4)=0.43607;
x(5)=1.33584;
x(6)=2.35060;

w(1)=0.00453;
w(2)=0.15706;
w(3)=0.724629;
w(4)=0.724629;
w(5)=0.157067;
w(6)=0.00453;

id_avg=0.021;
id_std=0.068; % percent

pcounter=1;
while PresDropSet<=PresDropMax
    error_f=0.1;
    f_iterations=0;
    while error_f>1e-4 && f_iterations<100
        error_f=0;
        for n=1:now
            id=id_avg*(1.4142*id_std*x(n)+1);

            ir=id/2; %cm
            TotalArea=2*pi*or*length*nof;
            A=TotalArea/nos;
            for m=1:2
                for i=1:noc
                    F(i)=z(i)*Ftotal(n,m);
                end

                for j=1:nos    %Do crossflow calculations
                    Ptotal(j)=0;
                    Rtotal(j)=0;
                    for i=1:noc   %Initialize with permeate vacuum
                        if j==1
                            P(j,i)=Q(i)*A*FeedPres*F(i)/Ftotal(n,m);
                            R(j,i)=F(i)-P(j,i);
                        else
                            P(j,i)=Q(i)*A*FeedPres*R(j-1,i)/Rtotal(j-1);
                            R(j,i)=R(j-1,i)-P(j,i);
                        end
                        Ptotal(j)=Ptotal(j)+P(j,i);
                        Rtotal(j)=Rtotal(j)+R(j,i);
                    end

                    c_error=0.1;
                    iterations=0;

                    while c_error>errormax && iterations<maxiterations   %Iterate crossflow solution
                        c_error=0;
                        for i=1:noc
                            if j==1
                                P_new=(Q(i)*A*FeedPres*F(i)/Rtotal(j))/(1+Q(i)*A*PermPres/Ptotal(j)+Q(i)*A*FeedPres/Ftotal(n,m));
                                R_new=F(i)-P_new;
                            else
                                P_new=(Q(i)*A*FeedPres*R(j-1,i)/Rtotal(j))/(1+Q(i)*A*PermPres/Ptotal(j)+Q(i)*A*FeedPres/Rtotal(j));
                                R_new=R(j-1,i)-P_new;
                            end
                            c_error=c_error+(abs(R_new-R(j,i))+abs(P_new-P(j,i)))/Ftotal(n,m);
                            iterations=iterations+1;
                            P(j,i)=P_new;
                            R(j,i)=R_new;
                        end

                        Ptotal(j)=0;
                        Rtotal(j)=0;
                        for i=1:noc
                            Ptotal(j)=Ptotal(j)+P(j,i);
                            Rtotal(j)=Rtotal(j)+R(j,i);
                        end
                    end
                end

                Sweeptotal=0;    %Sum permeates from last stage to first
                for i=1:noc
                    Sweep(i)=0;
                end

                for j=nos:-1:1
                    if j==nos
                        Ptotal(j)=Ptotal(j)+Sweeptotal;
                    else
                        Ptotal(j)=Ptotal(j)+Ptotal(j+1);
                    end

                    for i=1:noc
                        if j==nos
                            P(j,i)=P(j,i)+Sweep(i);
                        else
                            P(j,i)=P(j,i)+P(j+1,i); 
                        end
                    end
                end
                
                %Permeate flow 
             if m ==1 
               for j = 1:nos
                  for n = 1:now
                     for i=1:noc
                        Permf(n,i)=P(j,i);
                     end
                  end
               end
             end

                for j=1:nos
                    RetPres(j)=FeedPres;
                end

                d_maxiterations=500;
                c_error=0.1;   %Use direct-substitution method
                d_iterations=0;
                while c_error>errormax && d_iterations<=d_maxiterations
                    c_error=0;
                    functiontemp=Ftotal(n,m);
                    RetPres=RetPresDrop(R,Rtotal,RetPres);
                    for j=1:nos
                        Ptotal(j) = 0;
                        Rtotal(j) = 0;
                        P_sum(j) = 0;
                        for i=1:noc
                            for n = 1:now
                               P(j,i) = P_sum(j) + 1/sqrt(pi)*(w(n)*Permf(n,i)); %Summing up the flows from each fiber type(n)
                            end
                              Ptotal(j) = Ptotal(j) + P(j,i);
                              Rtotal(j) = Rtotal(j) + R(j,i);
                        end
                             for i=1:noc
                               y(j,i) = P(j,i)/Ptotal(j);
                             end
                          for i=1:noc
                            if j==1
                                   x(j,i) = F(i)/Rtotal(j);
                                   %y(j,i) = P(j,i)/Ptotal(j);
                                   J(j,i) = Q(i)*A*(RetPres(j)*x(j,i)-PermPres*y(j,i));
                                   R_new=  F(i) - J(j,i);
                                   P_new = P(j+1,i) + J(j,i);
%                                 P_new=(Q(i)*A*RetPres(j)*F(i)/Rtotal(j)+P(j+1,i)*(1+Q(i)*A*RetPres(j)/Rtotal(j)))/(1+Q(i)*A*PermPres/Ptotal(j)+Q(i)*A*RetPres(j)/Rtotal(j));
%                                 R_new=F(i)-(P_new-P(j+1,i));
                            elseif j==nos
                                   x(j,i) = R(j,i)/Rtotal(j);
                                  % y(j,i) = Sweep(i)/Ptotal(j);
                                   J(j,i) = Q(i)*A*(RetPres(j)*x(j,i)-PermPres*y(j,i));
                                   R_new=  R(j-1,i) - J(j,i);
                                   P_new = Sweep(i) + J(j,i);
%                                 P_new=(Q(i)*A*RetPres(j)*R(j-1,i)/Rtotal(j)+Sweep(i)*(1+Q(i)*A*RetPres(j)/Rtotal(j)))/(1+Q(i)*A*PermPres/Ptotal(j)+Q(i)*A*RetPres(j)/Rtotal(j));
%                                 R_new=R(j-1,i)-(P_new-Sweep(i));
                            else
                                   x(j,i) = R(j,i)/Rtotal(j);
                                   %y(j,i) = P(j,i)/Ptotal(j);
                                   J(j,i) = Q(i)*A*(RetPres(j)*x(j,i)-PermPres*y(j,i));
                                   R_new=  R(j-1,i) - J(j,i);
                                   P_new = P(j+1,i) + J(j,i);
%                                 P_new=(Q(i)*A*RetPres(j)*R(j-1,i)/Rtotal(j)+P(j+1,i)*(1+Q(i)*A*RetPres(j)/Rtotal(j)))/(1+Q(i)*A*PermPres/Ptotal(j)+Q(i)*A*RetPres(j)/Rtotal(j));
%                                 R_new=R(j-1,i)-(P_new-P(j+1,i));
                           end
                            c_error=c_error+(abs(R_new-R(j,i))+abs(P_new-P(j,i)))/Ftotal(n,m);
                            P(j,i)=P_new;
                            R(j,i)=R_new;
                          end

                        Ptotal(j)=0;
                        Rtotal(j)=0;
                        for i=1:noc
                            Ptotal(j)=Ptotal(j)+P(j,i);
                            Rtotal(j)=Rtotal(j)+R(j,i);
                        end
                        
                        if m==1
                            if j==nos
                                for i=1:noc
                                    Retf(n,i)=R(nos,i);
                                end
                                    Rtotaltemp=Rtotal(nos);
                            end
%                                 if j == 1
%                                   for i=1:noc
%                                     Permf(n,i)=P(j,i);
%                                   end
%                                     Ptotaltemp=Ptotal(1);
%                                 end              
                            if m==1
                                dimrec(n,j)=0;
%                                 dimperm(n,j) = 0;
                                for i=1:noc
                                    dimrec(n,j)=dimrec(n,j)+R(j,i);
%                                     dimperm(n,j) = dimperm(n,j) + P(j,i);
                                    dimdist(j)=j/nos;
                                end
                                dimxfrac(n,j)=R(j,1)/dimrec(n,j);
%                                 dimyfrac(n,j)=P(j,1)/dimperm(n,j);
                                dimrec(n,j)=dimrec(n,j)/Ftotal(n,1);
%                                 theta(n,j) = dimperm(n,j)/Ftotal(n,1);
                            end
                        end
                    end
                    d_iterations=d_iterations+1;
                end
                PresDrop(m)=(FeedPres-RetPres(nos))*14.7/76;

                Ftotal(n,2)=1.1*Ftotal(n,1);
            end

            if n==3
                xfrac_nv(pcounter)=Retf(3,1)/Rtotaltemp;
                recovery_nv(pcounter)=Rtotaltemp/Ftotal(3,1);
            end
            
            alpha=0.1*Ftotal(n,1)/(PresDrop(2)-PresDrop(1));
            error_f=error_f+abs(alpha*(PresDropSet-PresDrop(1))/Ftotal(n,1));
            Ftemp(n)=Ftotal(n,1);
            Ftotal(n,1)=Ftotal(n,1)+alpha*(PresDropSet-PresDrop(1));
            
        end
        f_iterations=f_iterations+1;
    end
    
    RetTotal=0;
   % FeedTotal=0;
    %PermTotal= 0;
    for i=1:noc
        Ret(i)=0;
        %Feed(i)=0;
        Perm(i)=0;
        %Permf_nos(i)=0;
        for n=1:now
            Ret(i)=Ret(i)+1/sqrt(pi)*(w(n)*Retf(n,i));
            %Perm(i)=Perm(i)+1/sqrt(pi)*(w(n)*Permf(n,i)); %% total permeate at stage 1
            %Feed(i)=Feed(i)+1/sqrt(pi)*(w(n)*Ftemp(n));
        end
        RetTotal=RetTotal+Ret(i);
        %FeedTotal=FeedTotal+Feed(i);
       % PermTotal=PermTotal + Perm(i);
    end
    
    FeedTotal=0;
    for n=1:now
         FeedTotal=FeedTotal+1/sqrt(pi)*(w(n)*Ftemp(n));
    end
    
   % MBCheck = - PermTotal + FeedTotal - RetTotal;
    xfrac_v(pcounter)=Ret(1)/RetTotal;
    recovery_v(pcounter)=RetTotal/FeedTotal;
    pcounter=pcounter+1;
    PresDropSet=PresDropSet*1.1;
     

end

function RetPres = RetPresDrop(R,Rtotal,RetPres)
global nos noc nof CompMW CompSigma FeedTemp FeedPres ...
        length ir functiontemp
for j=1:nos
    RetMu=0;

    for i=1:noc
        CompMu=((FeedTemp*CompMW(i))^0.5)/(CompSigma(i)^2);
        RetMu=RetMu+CompMu*R(j,i)/Rtotal(j);
    end

    RetMu=RetMu*2.6693e-6;
    if j==1
        deltaP=(8*RetMu*functiontemp*(76/FeedPres)*(FeedTemp/273.15)* ...
                  (length/nos))/(pi*(ir^4)*nof);
        deltaP=deltaP*7.5e-4;  %Pa to cmHg
        RetPres(j)=FeedPres-deltaP;
    else
        deltaP=(8*RetMu*Rtotal(j-1)*(76/RetPres(j-1))*(FeedTemp/273.15)* ...
                  (length/nos))/(pi*(ir^4)*nof);
        deltaP=deltaP*7.5e-4;  %Pa to cmHg
        RetPres(j)=RetPres(j-1)-deltaP;
    end
end
end

Answer 1

[Walter Roberson]

Around line 62 or so you have

for n=1:now

Within that for loop, around line 140 or so you have

                 if m ==1
                    for j = 1:nos
                        for n = 1:now
                            for i=1:noc
                                Permf(n,i)=P(j,i);
                            end
                        end
                    end
                end

Notice that you have the for n inside the for n . MATLAB does define what happens when you do this: the change to n will be in effect until the next outer for n starts, at which point n will be set to the "next" value of the outer loop as-if the change to n had not been done. But meanwhile until the bottom of the loop, n has changed. In particular at your line just inside the inner while loop you are having problems with,

                     functiontemp=Ftotal(n,m);

the n that is going to be used there is the n that was left-over from the inner for n loop, not the n associated with the outer for n loop.

Notice as well that the assignment to Permf is inside a triple for loop but is assigning to a variable using only two indices. The value that Permf(n,i) is going to be left with is going to be the last value assigned to it. You could achieve the same effect with

Permf(1:now, 1:noc) = P(nos, 1:noc)

and since you are not initializing Permf (growing it inside the loop) and it is always the same size, this in turn could be expressed as

Permf = repmat(P(nos, :), now, 1);

with no loop.

Are you sure that is what you want Permf to be assigned?

Any time you have a for loop that is not doing reduction (e.g., totaling the P values) and has more nesting levels than the number of indices of the output variable, then you should be optimizing to do only what would be left in the final result after the overly-nested assignments were complete.