[MATLAB数学相关] 偏微分方程的编程解法

[niedongbing]

function [T,r,t] = NDB(r0,rf,tf,it0,brf,N,I)
% r0,rf左右界限
% tf为时间下限
% it0为函数T的初值
% br0,brf分别为左边界值和右边界值
% N,I为r轴和时间轴区间数

dr = (rf-r0)/N; r = r0+[0:N]*dr;
dt = tf/I; t = [0:I]'*dt; %建立内节点

% I.C.初值
for n = 1:N+1
T(n,1) = it0(r(n));
end
% B.C.边界值
for i = 1:I+1
T( N+1,i) = brf(t(i));
end
% Newton Boundary (left)牛顿边值（采用中间差分牛顿边界得到）
m1 = dt/dr/dr; m2 = dr/dt; m3 = dt/dr/pi;
for i = 1:I
T(1,i+1) = (2*m1+1/m2)*T(2,i)+(1-2*m1-1/m2)*T(1,i)+ m3;
end

% 内节点的值（显示差分偏微分方程得到）
for i = 1:I
for n = 2:N
T(n,i+1) = (m1+1/(m2+(n-1)/m1))*T(n+1,i)+(1-2*m1-1/(m2+(n-1)/m1))*T(n,i)+m1*T(n-1,i);

end
end

以上为解方程编写的M函数
下面代入初边值调用函数(另一M函数)
function zp()
r0 = 1; rf = 1000; %左边界为1，右边界为1000
N = 999; I = 500; tf = 100; % r轴等分区间数为999，时间轴区间数为500
it0 = inline('0'); % 初始值

brf = inline('0'); %右边界值

[T,r,t] = NDB(r0,rf,tf,it0,brf,N,I);
% 以下是画图程序
figure
plot(r,T(r,1))
axis([0 100 0 0.1])
title('Solution')
xlabel('r')
ylabel('T')
hold on
plot(r,T(r,10))
hold on
plot(r,T(r,100))

上面的程序可以运行，但是结果和文献上不符合，而且上面程序中dr只能取1，取0.1或0.01是得不到图像的，请各位高手帮忙给看一看，多谢了！
方程如下：dT2/dr2+1/r*(dT/dr)=dT/dt 一维的轴坐标
初值：t=0, T=0
边界值：r=1,dT/dr=-1/(2*pi);
r=1000,T=0;
是不是自己差分的有问题，或是边界处理的不当？（若调用内置函数解此问题是完全可行的，已试过和文献吻合，但要进一步的结果却必须编程）
在差分r时，由于左边界为1，所以r=1+(n-1)*dr,是不是这个有问题了？
还请各位高手帮住解答，多谢！