用网格法计算透平叶片的温度分布的时候，一提高网格精度，右上角温度的分布就会不正常变化，再提高精度就不收敛了，想问一下大家怎么解决呢？
这是我的代码，截取叶片通风口右上角1\4的区域计算
m=input('please input the value of m='); %从上到下的行数
n=input('please input the value of n='); %从左到右的列数
To=1700;Ho=1000;Ti=400;Hi=250;
dx=3/(m-1);dy=3.5/(n-1);
lam=input('please input the value of lambda=');
t = zeros(m,n);
global value;
for r=(m*3/7+1):m
for c=1:(n*5/8-1) %冷却通道内部点的温度，不参与迭代。
t(r,c)=Ti;
end
end
for a=1:5000
v=max(t);
value=max(v);
s=t;
t(1,1)=lam*dy*(2*t(1,2))/(2*dx*(Ho*dx+lam*dy/dx+lam*dx/dy))+lam*dx*t(2,1)/(dy*(Ho*dx+lam*dy/dx+lam*dx/dy))+Ho*dx*To/(Ho*dx+lam*dy/dx+lam*dx/dy);
%左上外角点
for c=2:(n-1) %上部平直边界上的节点
t(1,c)=lam*dy*(t(1,c+1)+t(1,c-1))/(2*dx*(Ho*dx+lam*dy/dx+lam*dx/dy))+lam*dx*t(2,c)/(dy*(Ho*dx+lam*dy/dx+lam*dx/dy))+Ho*dx*To/(Ho*dx+lam*dy/dx+lam*dx/dy);
end
t(1,n)=lam*dy*(2*t(1,n-1))/(2*dx*(Ho*dx+lam*dy/dx+lam*dx/dy))+lam*dx*t(2,n)/(dy*(Ho*dx+lam*dy/dx+lam*dx/dy))+Ho*dx*To/(Ho*dx+lam*dy/dx+lam*dx/dy);
%右上外角点
for r=2:(m*3/7-1) %左边界上两角点之间的节点
t(r,1)=dy*(t(r-1,1)+t(r+1,1)+2*t(r,2))/(dx*(2*dy/dx+2*dx/dy));
end
for r=2:(m*3/7-1)
for c=2:(n-1) %一般节点part1
t(r,c)=dy*(t(r-1,c)+t(r+1,c)+t(r,c-1)+t(r,c+1))/(dx*(2*dy/dx+2*dx/dy));
end
end
t(m*3/7,1)=lam*dy*(2*t(m*3/7,2))/(2*dx*(Hi*dx+lam*dy/dx+lam*dx/dy))+lam*dx*t(m*3/7-1,1)/(dy*(Hi*dx+lam*dy/dx+lam*dx/dy))+Hi*dx*Ti/(Hi*dx+lam*dy/dx+lam*dx/dy);
%环通风口左上外角点
for c=2:(n*5/8-1) %环通风口横向平直边界上的节点
t(m*3/7,c)=lam*dy*(t(m*3/7,c-1)+t(m*3/7,c+1))/(2*dx*(Hi*dx+lam*dy/dx+lam*dx/dy))+lam*dx*t(m*3/7-1,c)/(dy*(Hi*dx+lam*dy/dx+lam*dx/dy))+Hi*dx*Ti/(Hi*dx+lam*dy/dx+lam*dx/dy);
end
t(m*3/7,n*5/8)=lam*dy*(2*t(m*3/7,n*5/8+1)+t(m*3/7,n*5/8-1))/(2*dx*(Hi*(dx+dy)/2+3*lam*dy/(2*dx)+3*lam*dx/(2*dy)))+lam*dx*(2*t(m*3/7-1,n*5/8)+t(m*3/7+1,n*5/8))/(2*dy*(Hi*(dx+dy)/2+3*lam*dy/(2*dx)+3*lam*dx/(2*dy)))+Hi*(dx+dy)*Ti/(2*(Hi*(dx+dy)/2+3*lam*dy/(2*dx)+3*lam*dx/(2*dy)));
%环通风口右上内角点
for r=(m*3/7+1):(m-1) %环通风口纵向平直边界上的节点
t(r,n*5/8)=lam*dx*(t(r-1,n*5/8)+t(r+1,n*5/8))/(2*dy*(Hi*dy+lam*dy/dx+lam*dx/dy))+lam*dy*t(r,n*5/8+1)/(dx*(Hi*dx+lam*dy/dx+lam*dx/dy))+Hi*dy*Ti/(Hi*dy+lam*dy/dx+lam*dx/dy);
end
for r=(m*3/7):(m-1)
for c=(n*5/8+1):(n-1) %一般节点part2
t(r,c)=dy*(t(r-1,c)+t(r+1,c)+t(r,c-1)+t(r,c+1))/(dx*(2*dy/dx+2*dx/dy));
end
end
t(m,n*5/8)=lam*dx*(2*t(m-1,n*5/8))/(2*dy*(Hi*dy+lam*dy/dx+lam*dx/dy))+lam*dy*t(m,n*5/8+1)/(dx*(Hi*dx+lam*dy/dx+lam*dx/dy))+Hi*dy*Ti/(Hi*dy+lam*dy/dx+lam*dx/dy);
%环通风口右下外角点
for c=(n*5/8+1):(n-1) %下边界上两角点之间的节点
t(m,c)=dy*(2*t(m-1,c)+t(m,c-1)+t(m,c+1))/(dx*(2*dy/dx+2*dx/dy));
end
for r=2:(m-1) %右边界上两角点之间的节点
t(r,n)=dy*(t(r-1,n)+t(r+1,n)+2*t(r,n-1))/(dx*(2*dy/dx+2*dx/dy));
end
t(m,n)=dy*(2*t(m-1,n)+2*t(m,n-1))/(dx*(2*dy/dx+2*dx/dy));
%右下外角点
M=max(t);
Max=max(M);
if abs(Max-value)<=1.0e-6
disp(['截面中的最高温度为：' num2str(value) 'K。']);
[row,col]=find(s==value);
disp(['最高温度点坐标为：(' num2str(row) ',' num2str(col) ')']);
disp('此时的温度分布为：');
disp(s);
break;
elseif a==5000
disp('方程不收敛。')
end
end