Решить задачу теплопроводности стального стержня

Введём регулярную сетку:
tn=nτ; xm=mh; umn=uxm, tn, n=0,…,N, m=0,…,M,
Получаем следующую разностную схему:
umn+1-umnτ=a2∙um-1n-2umn+um+1nh2
Тогда решение в следующий момент времени tn+1
umn+1=umn+a2τh2um-1n-2umn+um+1n
Условие устойчивости разностной схемы:
τ≤h22a2
Примем h=0.01, тогда пусть τ=h24a2=2.
M=lh;N=Tτ
При этом краевые условия:
u0n=0, n=0, 1, 2,…, N
uMn=T0, n=0, 1, 2,…,N
Начальные условия:
um0=T0xm2l,m=0, 1, 2,…,M
Основная программа (main2.m):
close all
clear all
clc
% Явная схема для уравнения теплопроводности на отрезке
l=2; T=60; T0=70;
lambda=45.4; c=460; ro=7200;
a2=lambda/(c*ro);

%Создание сетки:
h=0.01;% шаг по x
M=l/h+1;% число узлов по x
tau = round((h^2)/(4*a2));% шаг по t
N= T/tau+1;% число узлов по t:
x=zeros(M, 1);
for i=1:M x(i)=(i-1)*h; end
t=zeros(N,1);
for j=1:N t(j)=(j-1)*tau; end

% численное решение:
y= zeros(N,M);
% начальное условие:
for i=1:M
y(1,i)=T0*(x(i)^2)/l;
end
for j=1:(N-1)
% решение уравнения на внутренних узлах сетки:
for i=2:(M-1)
y(j+1,i)=y(j,i)+tau*a2*(y(j,i+1)-2*y(j,i)+y(j,i-1))/(h*h);
end
% граничные условия:
y(j+1,1) =0;
y(j+1,M) =T0;
end

figure
surf(x,t,y)
xlabel('x')
ylabel('t')
zlabel('y')
title('numerical solution')
Результаты выполнения: