Clear["Global`*"];

Problem

Solve numerically the differential equation

((∂u(x,t))/(∂t))=(((∂^2)u(x,t))/(∂(x^2)))

with the following initial condition

u(x,0)=(e^(-(x^2)))

and Dirichlet  boundary conditions

u(-∞,t)=0
u(+∞,t)=0

Exact solution

u0[x_,t_]=Exp[-x^2/(1+4t)]/Sqrt[1+4t];
sol=DSolve[{D[u[x,t],t]==D[u[x,t],x,x],u[x,0]==u0[x,0]},u,{x,t}]
Plot3D[u[x,t]/.sol,{x,-15,15},{t,0,10},PlotRange->All,PlotPoints->100]

Numerical solution using built-in function

{xmin, xmax}={-30, 30};
{tmin, tmax}={0,10};
numsol=NDSolve[{D[u[x,t],t]==D[u[x,t],x,x],u[x,0]==u0[x,0],u[xmin,t]==u0[xmin,0],u[xmax,t]==u0[xmax,0]},u,{x,xmin,xmax},{t,tmin,tmax}]
Plot3D[Evaluate[u[x,t]/.numsol],{x,xmin,xmax},{t,tmin,tmax},PlotRange->All]

Numerical solution using basic explicit methods

(* grid initialization *)
{nx,nt}={201,301};
dx=(xmax-xmin)/(nx-1);
dt=(tmax-tmin)/(nt-1);
Print["λ = ",λ=N[dt/dx^2]]
X = N[Range[xmin, xmax, dx]];
T = N[Range[tmin, tmax, dt]];

(* Initialization of the array with zeroes - Dirichlet's boundary conditions *)
v=ConstantArray[0.0,{nx,nt}];
error=ConstantArray[0.0,nt];

(* Initial state *)
Do[v[[i,1]]=u0[X[[i]],0],{i,2,nx-1}];
Do[v[[i,2]]=u0[X[[i]],0],{i,2,nx-1}];

method=1;
Which[
method==1,
Print["Euler explicit method - order 1, stable for λ = k/h^2 < 1/2"];
Do[v[[j,n+1]]=v[[j,n]]+λ(v[[j+1,n]]-2v[[j,n]]+v[[j-1,n]]),{n,1,nt-1},{j,2,nx-1}],
(* Leap-frog method - unstable *)
method==2,
Print["Leap-frog method - unstable"];
Do[v[[j,n+1]]=v[[j,n-1]]+2λ(v[[j+1,n]]-2v[[j,n]]+v[[j-1,n]]),{n,2,nt-1},{j,2,nx-1}]
];
Do[
error[[n]]=0.0;
Do[
error[[n]]=Max[error[[n]],Abs[v[[j,n]]- u0[X[[j]],T[[n]]]]],
{j,2,nx-1}
],
{n,1,nt-1}
];

(* fancy plotting *)
Manipulate[
n=Round[t/dt]+1;
ListLinePlot[Table[{X[[j]],v[[j,n]]},{j,1,nx}],
PlotRange->{-0.5,1.5}],
{t,tmin,tmax,dt}]

(* show maximal errors *)
ListLogPlot[{Table[{T[[n]],error[[n]]},{n,1,nt}]},
PlotRange->All]
Print["Maximal error: ",Max[error]]