computer code from R. M. Feenstra, “Path distributions for describing eigenstates of the harmonic oscillator and other 1-dimensional problems”,  arXiv:2306.11155 [quant-ph].

computations for Fig. 7, providing results of Eq. (36) using Eqs. (35), (34), (33) and (31)

in variable names below, an appended ‘p’ means prime (so Tp = T’ (T-prime), pp = p’ (p-prime), etc.)

turning points are when M (ω^2) (x^2)/2 = (n + 1/2) hbar ω, so x = Sqrt[(2n + 1) hbar / M ω]; perform xf integral over 5 of these ranges for either sign of xf

use specified nTp number of points for the integral over Tp (32 for fine computations, 2 for coarse)

compute results from p = 0 to 5 with specified delp spacing (0.1 for fine computations, 0.5 for coarse)

includes correction on first discrete point evaluated in pp integral, near the 1/Sqrt[(pp^2)-(xf^2)(ω^2)] divergence (integrate that term over the delpp interval)

Clear["Global`*"]; 
ω=1; T=N[32 π]; nTp=32;delTp=N[2π/nTp];
pmax=5; delp=0.1; np=Round[pmax/delp];
xfmax=N[5Sqrt[(2n+1)/ω]]; delxf=0.1; nxf=2Ceiling[xfmax/delxf]; delxfp=2 xfmax/nxf;
npp=Floor[Max[50,150Abs[xf]/(Sqrt[(2n+1)/ω])]];
prefact[ω_,T_]:=N[((((ω/(N[π]))))^(1/2))((Sqrt[Abs[Sin[T ω]]])/((2^n)Factorial[n](Sqrt[ 2 N[π] I ω ])))Exp[I ((n+1/2)ω T-π Floor[ω T/N[π]]/2)]];
x0[p_,xf_,ω_,T_]:=xf Cos[T ω]-Sign[p]Sin[T ω]√((((p/ω))^2)-(xf^2));
S[p_,xf_,ω_,T_]:=N[(Sin[T ω]/(2ω))( Cos[T ω]((p^2)-2(xf^2)(ω^2))+2 xf ω Sign[p]Sin[T ω]√((p^2)-(xf^2)(ω^2)))];
tab[n_,xfmax_,nxf_,delxfp_,np_]:=Table[{p=(ip+0.5)delp; p,
delxfp Total[Table[xf=-xfmax+(ixf-0.5) delxfp;Exp[-ω (xf^2)/2]HermiteH[n,(Sqrt[ω])xf](((2 N[π]))^(-1))delTp Total[Table[Tp=N[T+(iTp-0.5)delTp]; delpp=1/(Sqrt[Tp]npp);  
(Sqrt[Tp]/2)Re[prefact[ω,Tp]Total[Table[pp=p+(ipp+0.5)delpp; If[pp>ω Abs[xf]||pp<-ω Abs[xf],If[(pp>0 && (pp-ω Abs[xf])<=delpp)||(pp<0 && (-ω Abs[xf]-pp)<=delpp), Sqrt[(((Abs[pp]+0.5delpp))^2)-(xf^2)(ω^2)]HermiteH[n,√ω x0[pp,xf,ω,Tp]]Exp[-ω ((x0[pp,xf,ω,Tp])^2)/2]Exp[I S[pp,xf,ω,Tp]],delpp(Abs[pp]/Sqrt[(pp^2)-(xf^2)(ω^2)])HermiteH[n,√ω x0[pp,xf,ω,Tp]]Exp[-ω ((x0[pp,xf,ω,Tp])^2)/2]Exp[I S[pp,xf,ω,Tp]]],0],{ipp,-npp,npp-1}]]],{iTp,nTp}]],{ixf,nxf}]]},{ip,0,np-1}];
n=0;
tab0=tab[n,xfmax,nxf,delxfp,np];
fname="tab0.dat"; Print[fname];  stmp=OpenWrite[fname]; Write[stmp,tab0];  Close[stmp];
Print["total = ",2  delp Total[Table[tab0[[ip]][[2]],{ip,np}]]];
n=1;
tab1=tab[n,xfmax,nxf,delxfp,np];
fname="tab1.dat"; Print[fname];  stmp=OpenWrite[fname]; Write[stmp,tab1];  Close[stmp];
Print["total = ",2  delp Total[Table[tab1[[ip]][[2]],{ip,np}]]];
n=2;
tab2=tab[n,xfmax,nxf,delxfp,np];
fname="tab2.dat"; Print[fname];  stmp=OpenWrite[fname]; Write[stmp,tab2];  Close[stmp];
Print["total = ",2  delp Total[Table[tab2[[ip]][[2]],{ip,np}]]];
n=3;
tab3=tab[n,xfmax,nxf,delxfp,np];
fname="tab3.dat"; Print[fname];  stmp=OpenWrite[fname]; Write[stmp,tab3];  Close[stmp];
Print["total = ",2  delp Total[Table[tab3[[ip]][[2]],{ip,np}]]];

output for nTp=32; delp=0.1;

output for nTp=2; delp=0.5;