(*Griffiths Example 2.6. Run this cell.*)
ℏ:=1
m:=1
k:=1

(*Given equation for Ψ[x,0], and the result ϕ[k] found by using Plancherel's theorem. Run this cell.*)
Ψ0[x_]:=Piecewise[{{(1/(Sqrt[2a])),-a<x<a}},0]
ϕ[k_]:=(1/(Sqrt[π a]))((Sin[k a])/k)

Section 1 : Varying “a”

(*First, run this cell to see a plot of  Ψ[x,0] and ϕ[k] and for a given value of a =1.*)
a:=1
Plot[Ψ0[x],{x,-10,10},PlotRange->{-0.2,1},AxesLabel->{Style["x",Large],Style["Ψ[x,0]",Large]},Exclusions->None,AxesStyle-> Thick,TicksStyle->14]
Plot[ϕ[k],{k,-20,20},PlotRange->{-0.2,0.7},AxesLabel->{Style["k",Large],Style["ϕ[k]",Large]},AxesStyle-> Thick,TicksStyle->14]

Section 2 : Time - Dependence (We’ll actually use the PhET simulation for this part, but since I did code it up, here’s how you might use Mathematica to view the time dependence.)

(*An attempt to numerically approximate Ψ[x,t] without making Mathematica too angry. I've limited the range over which I'm integrating (not from - to + infinity in k) and instructed Mathematica to not do too many recursive calculations. *)
Ψ[x_,t_]:=(1/(Sqrt[2π]))NIntegrate[ϕ[k]Exp[-I(k x-ℏ (k^2)t/(2m))],{k,-30,+30},MaxRecursion->3]

Ψ[0,0]

(*Asking Mathematica to plot about 80 points will take significantly less time*)
DiscretePlot[{Abs[Ψ[x,0]]^2,Abs[Ψ[x,1]]^2},{x,-2a,+2a,0.05a},Filling->None,Joined->True,AxesLabel->{"x/a","Ψ^2[t]"}]

(*Below is the plot produced if I just ask Mathematica to "Plot" and set my NIntegrate to only bounds from -100 to +100 and MaxRecursion --> 3. It takes a long time. I have roughly approximated the original wavefunction!*)