clear;clc;

fcu=20;fc=9.6;ft=1.1;Ec=2.55*10^4;v=0.2;
h=800;b=100;s=100;Ast=615;Asvt=100;

G=Ec/(2*(1+v));
ets=0;et0=ft/Ec;
t=0;
while ets<=et0
      ets=ets+10^-6;
      t=t+1;
      phi=1/(3+1.8*b/h);
      phik=phi/(10*exp(1.3*b/h));
      thit(t)=20*ets/b*exp(1.3*b/h);
      T(t)=phik*h*b^2*G*thit(t);
end
Tcr=T(t);thcr=thit(t);

eds=10^-5;est=0;j(1)=0;
while eds<=0.0038&&est<=0.01
      eds=eds+0.0005;   
      t=t+1
      
      te=10;a=pi/4;zita=1;
      r=[0.4419*te,0.2648*te,0.1620*te,0,-0.1620*te,-0.2648*te,-0.4419*te];
      fda=0;mda=0;
      for i=1:7
          ed(i)=eds/te*(te-(te/2-r(i)));
          fd(i)=funfc(zita,ed(i));
          fda=fda+fd(i);
      end
      fdd=fda/7;
      for i=1:7
          md(i)=fd(i)/7*(te/2-r(i));
          mda=mda+md(i);
      end
      z=mda/fdd;
      edd=abs((te-z)/te*eds);
      ucor=2*(b+h)-8*z;
      Acor=(b-2*z)*(h-2*z);
      
      Tb=Acor*fdd*te*sin(2*a);
      evt=abs(Acor^2*fdd*eds/(ucor*Tb*tan(a))-edd);
      est=abs(Acor^2*fdd*eds/(ucor*Tb*cot(a))-edd);
      fsvt=abs(funfs(evt));fst=abs(funfs(est));
      
      te0=Ast*fst/(ucor*fdd)+Asvt*fsvt/(s*fdd);
      zita0=0.9/sqrt(1+400*(est+evt+edd));
      a0=acos(sqrt(Ast*fst/(ucor*fdd*te0)));
      
      while zita0<=zita
            zita=zita0;
            while a0<=a
                  a=a0;
                  
                  while abs(te-te0)>0.1
                        te=te0;
                        r=[0.4419*te,0.2648*te,0.1620*te,0,-0.1620*te,-0.2648*te,-0.4419*te];
                        fda=0;mda=0;
                        for i=1:7
                            ed(i)=eds/te*(te-(te/2-r(i)));
                            fd(i)=funfc(zita,ed(i));
                            fda=fda+fd(i);
                        end
                        fdd=fda/7;
                        for i=1:7
                            md(i)=fd(i)/7*(te/2-r(i));
                            mda=mda+md(i);
                        end
                        z=mda/fdd;
                        edd=(te-z)/te*eds;
                        ucor=2*(b+h)-8*z;
                        Acor=(b-2*z)*(h-2*z);
      
                        Tb=Acor*fdd*te*sin(2*a);
                        evt=Acor^2*fdd*eds/(ucor*Tb*tan(a))-edd;
                        est=Acor^2*fdd*eds/(ucor*Tb*cot(a))-edd;
                        fsvt=abs(funfs(evt));fst=abs(funfs(est));
                        
                        te0=Ast*fst/(ucor*fdd)+Asvt*fsvt/(s*fdd);
                        zita0=0.9/sqrt(1+400*(est+evt+edd));
                        a0=acos(sqrt(Ast*fst/(ucor*fdd*te0)));
                  end
            end
      end
      T(t)=Tb;
      thit(t)=eds/(te*sin(2*a));
end


plot(1000*thit,T/1000);
grid;
xlabel('θ单位：rad·m');
ylabel('T单位：kN*m');
title('T-θ关系图')
hold on
plot(1000*thcr,Tcr/1000,'r*')
text(1000*thcr,Tcr/1000,'Tcr')



function y=funfc(zita,ed)
fc=9.6;ft=1.1;Et=2.55*10^4;ec0=0.002;
et0=ft/Et;
et0=-et0;
if ed<=et0
   y=0;
end
if ed<=0&&ed>et0
   y=-Et*ed;
end
if 0<ed&&ed<=zita*ec0
   y=zita*fc*(2*ed/(zita*ec0)-(ed/(zita*ec0))^2);
end
if ed>zita*ec0&&ed<=0.0038
   y=zita*fc*(1-0.15*((ed-zita*ec0)/(0.0038-ec0)));
end
if ed>0.0038
    y=0;
end

function y=fun2(es)
Es=2*10^5;
if es<=-0.00105
    y=-210;
end
if es>-0.00105&&es<=0.00105
    y=Es*es;
end
if es>0.00105
    y=210;
end
y=abs(y)

这是我编制的一个土木工程上用的纯扭结构T-thita曲线，由于用了一些理论，需要同时调整三个参数，有办法让这个程序一定收敛吗？

或者谁有编好的扭转程序，拜托各位大神了！！