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

[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