Явна диференчна схема за едномерното уравнение на топлопроводността

heatEquation[l_,T_,h_,u0_,uL_,uR_,d_]:=(
τ=h^2/(3d);
n=Ceiling[l/h];
m=Ceiling[T/τ];
x=Table[(i-1)*h,{i,1,n+1}];
t=Table[(i-1)*τ,{i,1,m+1}];
y=Table[0,{n+1},{m+1}];

(*Initial condition*)
For[i=1,i≤n+1,i++,
y[[i,1]]=u0[x[[i]]];
];

(*Iterate over time*)
For[j=1,j≤m, j++,
	(*Boundary conditions*)
	y[[1,j+1]]=uL[t[[j]]];
	y[[n+1,j+1]]=uR[t[[j]]];

	(*Internal values from the main PDE*)
	For[i=2, i≤n, i++,
		y[[i,j+1]]=d τ/h^2y[[i-1,j]]+(1-2d τ/h^2)*y[[i,j]]+d τ/h^2y[[i+1,j]]
	];
];

{x,t,y}
)

Тестване на метода върху задачата (Subscript[u, t])-(Subscript[u, xx])=0, 0<x<1, 0<t≤0.1, u(x,0)=sin(2πx), u(0,t)=u(1,t)=0.

u0[x_]:=Sin[2Pi*x]
uL[x_]:=0
uR[x_]:=0
exactSol[x_,t_]:=E^(-4Pi^2*t)*Sin[2Pi*x]
apprSol = heatEquation[1,0.5,0.01,u0,uL,uR,1];
apprX = apprSol[[1]];
apprT = apprSol[[2]];
apprY = apprSol[[3]];

Manipulate[Show[Plot[exactSol[x, apprT[[j]]],{x,0,1},PlotStyle->Red, PlotRange->{{0,1},{-1,1}}],ListPlot[Transpose[{apprX,apprY[[All,j]]}]]],
 {j,1,Length[apprT]}]

scheme=(Blend[{RGBColor[1,0,0],RGBColor[0.75,0,0],RGBColor[1,0.55,0],White,RGBColor[0,0.55,1],RGBColor[0,0,0.75], RGBColor[0,0,1]},Rescale[#1,{-1,1}]]&);
Animate[
  ListDensityPlot[Flatten[{Transpose[{apprX,Table[0,{Length[apprX]}], apprY[[All,j]]}], Transpose[{apprX,Table[0.1,{Length[apprX]}], apprY[[All,j]]}]},1],
AspectRatio->Automatic,
ColorFunction->scheme,
ColorFunctionScaling->False,

PlotLegends->Placed[BarLegend[{scheme,{-1,1}},LabelStyle->Directive[Black,FontSize->12],LegendLabel->Placed["Temperature", Left],LegendMarkerSize->150],{Bottom,Center}]],
 {j,1,Length[apprT],10}]