有大侠可以帮忙一下吗，为什么我的编程一直在循环，是不是while出问题了，求帮助啊

chinasteel

b=1000
h=2000
fc=14.3
ec=0.00164
ft=1.43
et=65*10^(-6)*ft^0.54   
at=0.312*ft^2   
n=50     
fy=360
Es=2*10^5
es=fy/Es

  r=zeros(1,100)  
for k=1:1:100
     o=(5e-4)*k
     x1=0,x2=h,df=10
  while abs(df)>10
    x0=(x1+x2)/2,x=h-x0  

  s=zeros(1,(round (x0/n)-1))  ,t=zeros(1,(round (x0/n)-1))
  for p=1:1:(round (x0/n)-1)
    i=n*p
    e1=-((i-n/2)-x0)*o        
    e2=-(i-x0)*o               
      if i==50 | i==1950
            steelarea=8844
      elseif i==150 | i==250 | i==350 | i==450 | i==550 | i==650 | i==750 | i==850 | i==950 | i==1050 | i==1150 | i==1250 | i==1350 | i==1450 | i==1550 | i==1650 | i==1750 | i==1850
            steelarea=1608      
      else
            steelarea=0
      end

         if e1<=ec
             concretecompression=fc*(0.03*(e1/ec)^3-1.06*(e1/ec)^2+2.03*(e1/ec))      
         else
             concretecompression=(e1/ec)*fc/(1.36*(e1/ec-1)^2+e1/ec)
         end
         
         if e2<=es                                       
            steelcompression=Es*e2
         elseif e2<0.01
            steelcompression=fy
         else
            steelcompression=0
         end         
   
      
         df1=(concretecompression*b*n+steelcompression*steelarea)/1000
         s(p)=df1
         dM1=-df1*(i-x0)/1000
         t(p)=dM1
   end
      
  

   y=zeros(1,(round (x/n)-1))  ,z=zeros(1,(round (x/n)-1))
   for q=1:1:(round (x/n)-1)
     j=n*q
     e3=((j-n/2)-x)*o
     e4=(j-x)*o

      if j==50 | j==1950
            steelarea=8844
      elseif j==150 | j==250 | j==350 | j==450 | j==550 | j==650 | j==750 | j==850 | j==950 | j==1050 | j==1150 | j==1250 | j==1350 | j==1450 | j==1550 | j==1650 | j==1750 | j==1850
            steelarea=1608      
      else
            steelarea=0
      end

       if e3<=et
          concretetension=(1.2*(e3/et)-0.2*(e3/et)^6)*ft
       else
          concretetension=(e3/et)*ft/(at*(e3/et-1)^1.7+e3/et)
       end
  
      if e4<=es
          steeltension=Es*e4
      elseif e4<0.01
          steeltension=fy
      else
            steeltension=0
      end

         df2=(concretetension*b*n+steeltension*steelarea)/1000   
         y(q)=df2
         dM2=-df2*(j-x)/1000
         z(q)=dM2
  end
         f1=sum(r(:))
         M1=sum(s(:))
         f2=sum(y(:))
         M2=sum(z(:))         
        df=f1-f2
        M=M1+M2


    if df>0
        x1=x0
     else
        x2=x0
    end   
end
r(k)=M  
end
r
o=5e-4:5e-4:5e-2
plot(0,r)