我现在编写了一个求解一维热传导的偏微分方程，调程序都调了好多天了 不知道问题在哪儿，求各位高手帮忙看一下好么？ 我在此表示万分感谢

需求解的方程看图片，我的程序如下

%-------------------1.主函数----------------------

function heat_conductionZ()

global M;

global Fo;

global T;

%初始化数据

c=1.48;

p=1714.8;

s=0.535

options={'空间半径ro','空间网格数N' ,'时间网格数M','对流传热系数h','稳定条件的值Fo(取值必须小于0.5)',};

topic='seting';

lines=1;%采用plot绘线色

def={'0.032','10','100','61.27','0.5'};

b=inputdlg(options,topic,lines,def);%创建输入对话

ro=eval(b{1});%eval执行字符串

N=eval(b{2});

M=eval(b{3});

h=eval(b{4});

Fo=eval(b{5});

%***************************************************

a=s/(p*c)

dr=ro/N            %空间步长

r=0:dr:ro

r=r'

ot=Fo*dr^2/a       %时间步长

tm=M*ot            %热传导的总时间tm

t=0:ot:tm

t=t'

%计算初值和边值

T=zeros(N+1,M+1)

Ti=init_fun(r)                  %fun为目标函数的表达式字符串或MATLAB自定义函数的函数柄

To=border_funo(t);

Te=border_fune(t);

T(:,1)=Ti                       %t=0时的初始温度

T(1,:)=To                       %r=0时的温度值

T(N+1,:)=Te                     %r=r0时的温度

%用差分法求出温度T与半径r、时间t的关系

for k=1:M

m=2

while m<=N&&Fo<=(ro/(2*r+dr))

T(m,k+1)=(Fo+(dr*Fo)/ro)*(T(m+1,k)+(1-(dr*Fo)/r+2*Fo)*T(m,k))+Fo*T(m-1,k)

m=m+1

end;

end;

%设置立体网格

for i=1:M+1

X(:,i)=r

end;

for j=1:N+1

Y(j,:)=t

end

mesh(X,Y,T);

view([1 -1 1]);

xlabel('r');

ylabel('t');

zlabel('T');

%-------------------2.子函数----------------------

function y=border_funo(t)

global M;

global Fo;

global T;

m=1;

k=1;

while(k<=M)

T(m,k+1)=(1-2*Fo)*T(m,k)+2*Fo*T(m+1,k);

k=k+1;

y

end;

return

%------------------------------------------------------

%-------------------3.子函数----------------------

function y=border_fune(t)

global M;

global Fo;

global T;

global N;

Tf=40;

A=(h*ot)/(p*c*dr);

for k=1:M

m=N;

while (1-A-Fo)>=0

T(m+1,k+1)=A*Tf+(1-A-Fo)*T(m+1,k)+Fo*T(m,k);

y=T(m+1,k);

end;

end;

return

%------------------------------------------------------

%------------------4.子函数--------------------------

function y=init_fun(r)%初值条件t=0

y=14;

return

2014-3-8 19:58 上传