Clear["Global`*"];

Problem

Solve numerically the differential equation (in atomic units ℏ = 1, (Subscript[m, e]) = 1)

i((∂ψ(x,t))/(∂t))=-(1/(2μ))(((∂^2)ψ(x,t))/(∂(x^2)))+V(x)ψ(x,t)

with the following initial condition

ψ(x,0)=(1/(((2π (σ^2)))^(1/4)))(e^(-(((x-(Subscript[x, 0])))^2)/(4(σ^2)) + i p x))

and Dirichlet  boundary conditions

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

Exact solution for V(x) = 0

Initial normalized Gaussian packet :

u0[x_,x0_,σ_,p_]=1/(2π σ^2)^(1/4)Exp[(-(x-x0)^2/(4σ^2)+I p x)];
Assuming[σ>0,Integrate[u0[x,x0,σ,p]*u0[x,x0,σ,-p],{x,-Infinity,Infinity}]]

sol=DSolve[{I D[u[x,t],t]==-D[u[x,t],x,x]/(2μ),u[x,0]==u0[x,x0,σ,p]},u,{x,t}]
pini=3;mass=1;xini=-10;σini=1;
xmin=-15; xmax=15; tmax=5;
Plot3D[Abs[u[x,t]/.sol/.{σ->σini,x0->xini,p->pini,μ->mass}]^2,{x,xmin,xmax},{t,0,tmax},PlotRange->All,PlotPoints->100]
{Plot3D[Re[u[x,t]/.sol/.{σ->σini,x0->xini,p->pini,μ->mass}],{x,xmin,xmax},{t,0,tmax},PlotRange->All,PlotPoints->100],
Plot3D[Im[u[x,t]/.sol/.{σ->σini,x0->xini,p->pini,μ->mass}],{x,xmin,xmax},{t,0,tmax},PlotRange->All,PlotPoints->100]}

Numerical solution using built-in function

{xmin, xmax}={-20, 20};
{tmin, tmax}={0,5};
numsol=NDSolve[{I D[u[x,t],t]==-D[u[x,t],x,x]/(2mass),u[x,tmin]==u0[x,xini,σini,pini],u[xmin,t]==0,u[xmax,t]==0},u,{x,xmin,xmax},{t,tmin,tmax}]
Plot3D[Evaluate[Abs[u[x,t]/.numsol]^2],{x,xmin,xmax},{t,tmin,tmax},PlotRange->All,PlotPoints->100]

Numerical solution using Euler explicit method, unstable

(* grid initialization *)
{nx,nt}={201,501};
dx=(xmax-xmin)/(nx-1);
dt=(tmax-tmin)/(nt-1);
Print["λ/2m = ",λ=N[dt/dx^2/(2mass)]]
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]],xini,σini,pini],{i,2,nx-1}];

Print["Euler explicit method - order 1, unstable"];
Do[v[[j,n+1]]=v[[j,n]]+λ I(v[[j+1,n]]-2v[[j,n]]+v[[j-1,n]]),{n,1,nt-1},{j,2,nx-1}]

Do[
error[[n]]=0.0;
Do[
error[[n]]=Max[error[[n]],Abs[v[[j,n]]- u[x,t]/.sol/.{σ->σini,x0->xini,p->pini,μ->mass,x->X[[j]],t->T[[n]]}]],
{j,2,nx-1}
],
{n,1,nt-1}
];

(* fancy plotting *)
Manipulate[
n=Round[t/dt]+1;
ListLinePlot[Table[{X[[j]],Abs[v[[j,n]]^2]},{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]]