h=0.01;
τ=0.01;
l=1;
T=0.1;
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]]=Sin[2 Pi  (i-1)h]
];

(*Generate the matrix of the linear system. This can be done outside the for-loop, since the
matrix is constant*)
A=Table[0,{n+1},{n+1}];
A[[1,1]]=1;
A[[n+1,n+1]]=1;
For[i=2,i<n+1,i++,
A[[i,i-1]]=-1/(2h^2);
A[[i,i]]=1/τ+1/h^2;
A[[i,i+1]]=-1/(2h^2);
];

(*Iterate over the time nodes*)
For[j=1,j≤m,j++,
(*Generate the right-hand side vector of the linear system*)
b=Table[0,{n+1}];
For[k=2,k≤n,k++,
b[[k]]=((y[[k,j]])/τ)+((y[[k-1,j]]-2y[[k,j]]+y[[k+1,j]])/(2h^2))
];
yNew=LinearSolve[A,b];
y[[All,j+1]]=yNew
]

Manipulate[
ListPlot[Table[{x[[i]],y[[i,j]]},{i,1,Length[x]}],PlotRange->{{0,1},{-1,1}}],
{j,1,Length[t],1}]