heatEquation[l_,T_,h_,f_,u0_,uL_,uR_,d_]:=(
τ=((h^2)/(3d));
n=Ceiling[l/h];
m=Ceiling[T/τ];
t=Table[(j-1)τ,{j,1,m+1}];
x=Table[(i-1)h,{i,1,n+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 the time nodes*)
For[j=1,j≤m,j++,
(*Compute the internal values from the main PDE*)
For[i=2,i≤n,i++,
y[[i,j+1]]=(1-2((d τ)/(h^2)))y[[i,j]]+((d τ)/(h^2))(y[[i-1,j]]+y[[i+1,j]])+τ f[x[[i]],t[[j]]]
];
(*Compute the boundary values from the boundary conditions*)
y[[1,j+1]]=uL[t[[j+1]]];
y[[n+1,j+1]]=uR[t[[j+1]]]
];
{x,t,y}
)

f[x_,t_]:=0
u0[x_]:=Sin[2Pi x]
uL[t_]:=0
uR[t_]:=0
resX=heatEquation[1,0.1,0.01,f,u0,uL,uR,1][[1]];
resT=heatEquation[1,0.1,0.01,f,u0,uL,uR,1][[2]];
resY=heatEquation[1,0.1,0.01,f,u0,uL,uR,1][[3]];
Manipulate[
Show[
ListPlot[Table[{resX[[i]],resY[[i,j]]},{i,1,Length[resX]}],PlotRange->{{0,1},{-1,1}}],
Plot[(E^(-4Pi^2resT[[j]]))Sin[2Pi x],{x,0,1},PlotStyle->Red]
],
{j,1,Length[resT],1}]

Manipulate[
Show[
ListPlot[Table[{resX[[i]],Abs[resY[[i,j]]-(E^(-4Pi^2resT[[j]]))Sin[2Pi x[[i]]]]},{i,1,Length[resX]}]]
],
{j,1,Length[resT],1}]