[求助]加入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

[chytemp]

我运行的结果除了可以得出T=1073时计算的数值，还有如下信息，但看不大懂提示哪里错误，是否T2=solve(f)命令有问题？

Warning: Out of range or non-integer values truncated during conversion to character.
(Type "warning off MATLAB:nonIntegerTruncatedInConversionToChar" to suppress this warning.)
> In C:\MATLAB6p5\toolbox\symbolic\solve.m at line 80
In C:\MATLAB6p5\work\trychy.m at line 57
??? Error using ==> sym/sym (char2sym)
{

Error in ==> C:\MATLAB6p5\toolbox\symbolic\@sym\sym.m
On line 92 ==> S = char2sym(x);

Error in ==> C:\MATLAB6p5\toolbox\symbolic\solve.m
On line 97 ==> vars = ['[' findsym(sym(eqns),neqns) ']'];

Error in ==> C:\MATLAB6p5\work\trychy.m
On line 57 ==> T2=solve(f);

[wxw0827]

T =
1073
T2 =
672.2851
T4 =
372.5158
T1 =
560.2376
T3 =
447.0190
N =
0.6736
P =
4.3339e+003
Warning: List of equations is empty.
> In solve at 87
In trychy at 61
??? Error using ==> maple
Error, (in linsolve) 1st argument fails to evaluate to a matrix

Error in ==> sym.mldivide at 28
X = maple('linsolve',char(A),char(B),'''_rank''');

Error in ==> sym.mrdivide at 27
X = (A.'\B.').';

Error in ==> trychy at 62
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-
>>

[chytemp]

是只能得到T=1073的结果，其他的循环数值没有，后来问题算是初步解决了，用fzero代替fsolve语句解方程，虽然没用循环但是单次运行还是可以的，现在又有一个新的问题很棘手。
我要做一个三维的网格图（mesh)，来看N=N(m3, m4)随m3和m4两个变量的变化情况。程序我贴在楼下，在command window可以看出，当输入ezmesh(@drawmesh,[0.006,0.007,0.01,0.02])命令时，程序是运行的（虽然很慢很慢，在我laptop上大概要等一个小时才运行完），计算结果是有的，但是图却出不来，而且有如下错误提示：

??? Error using ==> sym.maple
Error, (in evalf/ln) integer too large in context

Error in ==> sym.double at 23
X = reshape(str2num(map2mat(char(maple('evalf',S(,d)))),siz);

Error in ==> lambertw at 15
W = double(lambertw(sym(k)));

Error in ==> drawmesh>@(T2)(-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)-(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))))) at 50
>>f=@(T2)(-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*
Error in ==> fzero at 497
fb = FunFcn(b,varargin{:});

Error in ==> drawmesh at 51
T2=fzero(f,1000)

Error in ==> specgraph\private\ezplotfeval>applyfun at 85
z(i) = feval(f,x(i),y(i));

Error in ==> specgraph\private\ezplotfeval at 70
z = applyfun(x,y);

Error in ==> ezgraph3>ezeval at 619
X = ezplotfeval(fx,arg1,arg2);

Error in ==> ezgraph3>surfplot at 546
u = ezeval(F,var,X,Y);

Error in ==> ezgraph3 at 49
[dummy,h] = surfplot(f,domain,surfstyle,cax,Npts,fixdomain,flabel,fargs);

Error in ==> ezmesh at 66
h = ezgraph3('mesh',args{:});

不知道是什么原因? 请各位指点一下，万分感谢！

[chytemp]

程序虽然有些长，我想直接贴出来比较方便吧，如下：

function N=drawmesh(m4,m5)

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=14500;
T=1073;
%T=1173, i=17000
%T=1273, i=18500
% i can be 2000A/m^2 or 14500A/m^2 (where the power density of the fuel cell attains its maximum)

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;

x=1.2;

f=@(T2)(-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)-(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)))));
T2=fzero(f,1000)

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)

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

% to plot N-i-T type ezmesh(@drawmesh,[14500,17000,1073,1273]) in the command window