[求助]加入for循环后为何计算出错？

[chytemp]

编了个程序，其中变量T如手动设定一个数值，可正常计算，但若要实现变量T每隔10递增而采用了含for简单循环的语句（for T=1073:10:1173）时，程序运行时就没动静了，只能运行T=1073的数值，不知道是不是因为方程太大（是个非线性方程，没有解析解，只有数值解）没办法循环运算？
有请各位帮忙指点一下，万分感谢！

程序如下：
function T2=trychy(T)
syms n R F T0 ra rc Ea Ec ila ilc i p1 p2 p3 p0 Eel o L s h k T x m3 m4 m5 T2 T41 T42

n=2;
R=8.314;
F=96485;
T0=298;
ra=5.5*10^8;
rc=7.0*10^8;
Ea=1.0^10^5;
Ec=1.2*10^5;
ila=2.99*10^4;
ilc=2.16*10^4;
p1=0.97;
p2=0.21;
p3=0.03;
p0=1;
Eel=8*10^4;
o=3.6*10^7;
L=20*10^(-6);
s=-55.5666;
h=-248303;
k=1/100;

i=5000;

for T=1073:10:1173

g=h-T*s;

i0a=ra*(p1/p0)*(p3/p0)*exp(-Ea/(R*T));
i0c=rc*(p2/p0)^(1/4)*exp(-Ec/(R*T));

E0=-g/(n*F);
E=E0+R*T/(n*F)*log(p1*p2^(1/2)/p3);


d1=2*n*asinh(i/(2*i0a))+2*n*asinh(i/(2*i0c))-log(1-i/ila)-log(1-i/ilc)+i*n*F*L/(o*R)*exp(Eel/(R*T));
nint=R*T*d1/(n*F);

Pc=i*E-i*R*T*d1/(n*F)-(i*n*F*k)/(R*T*d1)*(E-R*T*d1/(n*F))^2;
Nc=Pc/(-i*h/(n*F));

V=E-R*T*d1/(n*F);

m3=1.1;
m4=0.05;
m5=0.01;

x=1.2;

T41=-T0*(exp(-(T2*x-T2-log((T*x-T2)/x/(T-T2))*x*i*m4+log((T*x-T2)/x/(T-T2))*x*i*m4*Nc+log((T*x-T2)/x/(T-T2))*x*m5*T-log((T*x-T2)/x/(T-T2))*x*m5*T0)/m3/x/(-i*m4+i*m4*Nc+m5*T-m5*T0))-1)/(-exp(-(T2*x-T2-log((T*x-T2)/x/(T-T2))*x*i*m4+log((T*x-T2)/x/(T-T2))*x*i*m4*Nc+log((T*x-T2)/x/(T-T2))*x*m5*T-log((T*x-T2)/x/(T-T2))*x*m5*T0)/m3/x/(-i*m4+i*m4*Nc+m5*T-m5*T0))+x);
T42=T0*(lambertw(exp((-m3*T*T2*x^2+m3*T^2*x^2+m3*T2^2*x-log((T*x-T2)/x/(T-T2))*T*T2*x^2-T*T2*x^2+log((T*x-T2)/x/(T-T2))*T^2*x^2-log((T*x-T2)/x/(T-T2))*T2^2-T*T2+2*T*T2*x-log((T*x-T2)/x/(T-T2))*T^2*x+log((T*x-T2)/x/(T-T2))*T*T2-m3*T*T2-2*m3*T*T2*x+m3*T^2*x+m3*T2^2+log((T*x-T2)/x/(T-T2))*T2^2*x)/m3/(T^2*x^2-T^2*x-T*T2*x^2+T*T2+T2^2*x-T2^2))/(x-1))*x-lambertw(exp((-m3*T*T2*x^2+m3*T^2*x^2+m3*T2^2*x-log((T*x-T2)/x/(T-T2))*T*T2*x^2-T*T2*x^2+log((T*x-T2)/x/(T-T2))*T^2*x^2-log((T*x-T2)/x/(T-T2))*T2^2-T*T2+2*T*T2*x-log((T*x-T2)/x/(T-T2))*T^2*x+log((T*x-T2)/x/(T-T2))*T*T2-m3*T*T2-2*m3*T*T2*x+m3*T^2*x+m3*T2^2+log((T*x-T2)/x/(T-T2))*T2^2*x)/m3/(T^2*x^2-T^2*x-T*T2*x^2+T*T2+T2^2*x-T2^2))/(x-1))-1)/(lambertw(exp((-m3*T*T2*x^2+m3*T^2*x^2+m3*T2^2*x-log((T*x-T2)/x/(T-T2))*T*T2*x^2-T*T2*x^2+log((T*x-T2)/x/(T-T2))*T^2*x^2-log((T*x-T2)/x/(T-T2))*T2^2-T*T2+2*T*T2*x-log((T*x-T2)/x/(T-T2))*T^2*x+log((T*x-T2)/x/(T-T2))*T*T2-m3*T*T2-2*m3*T*T2*x+m3*T^2*x+m3*T2^2+log((T*x-T2)/x/(T-T2))*T2^2*x)/m3/(T^2*x^2-T^2*x-T*T2*x^2+T*T2+T2^2*x-T2^2))/(x-1))*x-x-lambertw(exp((-m3*T*T2*x^2+m3*T^2*x^2+m3*T2^2*x-log((T*x-T2)/x/(T-T2))*T*T2*x^2-T*T2*x^2+log((T*x-T2)/x/(T-T2))*T^2*x^2-log((T*x-T2)/x/(T-T2))*T2^2-T*T2+2*T*T2*x-log((T*x-T2)/x/(T-T2))*T^2*x+log((T*x-T2)/x/(T-T2))*T*T2-m3*T*T2-2*m3*T*T2*x+m3*T^2*x+m3*T2^2+log((T*x-T2)/x/(T-T2))*T2^2*x)/m3/(T^2*x^2-T^2*x-T*T2*x^2+T*T2+T2^2*x-T2^2))/(x-1)));


f=subs(T41-T42);
T2=solve(f);
T4=-T0*(exp(-(T2*x-T2-log((T*x-T2)/x/(T-T2))*x*i*m4+log((T*x-T2)/x/(T-T2))*x*i*m4*Nc+log((T*x-T2)/x/(T-T2))*x*m5*T-log((T*x-T2)/x/(T-T2))*x*m5*T0)/m3/x/(-i*m4+i*m4*Nc+m5*T-m5*T0))-1)/(-exp(-(T2*x-T2-log((T*x-T2)/x/(T-T2))*x*i*m4+log((T*x-T2)/x/(T-T2))*x*i*m4*Nc+log((T*x-T2)/x/(T-T2))*x*m5*T-log((T*x-T2)/x/(T-T2))*x*m5*T0)/m3/x/(-i*m4+i*m4*Nc+m5*T-m5*T0))+x);


T=double(T)
T2=double(T2)
T4=double(T4)


T1=T2/x
T3=x*T4

u=(T-T1)/(T-T2);
v=(T3-T0)/(T4-T0);


Nt=1-T4*(x-1)/(m3*log(v))*((i*m4*(1-Nc)-m5*(T-T0))^(-1)-log(u)/(T2*(1-1/x)));
Pt=(-i*h/(n*F)*(1-Nc)+h*m5/(n*F*m4)*(T-T0))*Nt;

N=Nc+Nt*(1-Nc-m5*(T-T0)/(m4*i))
P=(-i*h/(n*F))*N

end