Note new version of Mathematica knows how to Integrate AiryAi[x]^2,
so you can change the below NIntegrates to Integrates; however is slower
---
a[1]=-2.3
a[2]=-4.1
a[3]=-5.5
a[4]=-6.8
a[5]=-7.9
a[6]=-9
a[7]=-10
a[8]=-11
a[16]=-17.6
a[32]=-28.2
For[i=1,i<9,i++,a[i]=x /. N[FindRoot[AiryAi[x],{x,a[i]},AccuracyGoal->13],12] ]
a[16]=x /. N[FindRoot[AiryAi[x],{x,a[16]},AccuracyGoal->13],12]
a[32]=x /. N[FindRoot[AiryAi[x],{x,a[32]},AccuracyGoal->13],12]
aint[1]=NIntegrate[AiryAi[x]^2,{x,a[1],5}]
n[1]=1/Sqrt[aint[1]]
For[i=2,i<9,i++,aint[i]=aint[i-1]+NIntegrate[AiryAi[x]^2,{x,a[i],a[i-1]}] ]
For[i=2,i<9,i++,n[i]=1/Sqrt[aint[i]] ]
aint[16]=aint[8]+NIntegrate[AiryAi[x]^2,{x,a[16],a[8]}]
n[16]=1/Sqrt[aint[16]]
aint[32]=aint[8]+NIntegrate[AiryAi[x]^2,{x,a[32],a[16]}]
n[32]=1/Sqrt[aint[32]]

Plot[{(n[1]*AiryAi[x+a[1]])^2-a[1],(n[2]*AiryAi[x+a[2]])^2-a[2],(n[4]*AiryAi[x+a[4]])^2-a[4],(n[8]*AiryAi[x+a[8]])^2-a[8],x,
        -a[1],-a[2],-a[4],-a[8]},{x,0,15}, \
PlotRange->{0,12},PlotStyle->{{RGBColor[0,0,0]},{RGBColor[0,0,0]},
        {RGBColor[0,0,0]},{RGBColor[0,0,0]},{RGBColor[0,0,1]},{RGBColor[1,0,0]},     
         {RGBColor[1,0,0]},{RGBColor[1,0,0]},{RGBColor[1,0,0]}}]
Plot[(n[16]*AiryAi[x+a[16]])^2,{x,0,21}]
Plot[(n[32]*AiryAi[x+a[32]])^2,{x,0,31},PlotRange->{0,.22},AspectRation->.4]
Show[%,Plot[.095/Sqrt[-a[32]-x],{x,0,27.5},PlotStyle->{RGBColor[1,0,0]}] ]
For[i=1,i<9,i++,f[i,x_]=n[i]*AiryAi[x+a[i]] ]
f[16,x_]=n[16]*AiryAi[x+a[16]]
f[32,x_]=n[32]*AiryAi[x+a[32]]
-----
a[8] is n=287
---
a[1]=-120.2
a[2]=-120.5
a[3]=-120.8
a[4]=-121.1
a[5]=-121.4
a[6]=-121.65
a[7]=-121.95
a[8]=-122.2
a[9]=-122.5
a[10]=-122.8
a[11]=-123.1
a[12]=-123.35
a[13]=-123.65
a[14]=-123.9
a[15]=-124.2
For[i=1,i<16,i++,a[i]=x /. N[FindRoot[AiryAi[x],{x,a[i]},AccuracyGoal->13],12] ]
temp[x_]=Sum[AiryAi[x+i]^2,{i,-120,0,5}]
aint[1]=NIntegrate[temp[x],{x,0,5}]+NIntegrate[AiryAi[x]^2,{x,a[1],-120}]
n[1]=1/Sqrt[aint[1]]
For[i=2,i<16,i++,aint[i]=aint[i-1]+NIntegrate[AiryAi[x]^2,{x,a[i],a[i-1]}];
n[i]=1/Sqrt[aint[i]]]
n0=1/Sqrt[NIntegrate[AiryAi[x]^2,{x,-2.33810,6}]]
For[i=1,i<16,i++,b[i]=n[i]*n0*NIntegrate[AiryAi[x]*AiryAi[x+a[i]-a[8]],{x,-2.33810,6}]]
bs=0
For[i=1,i<16,i++, Print[{i,b[i]}]; bs=bs+b[i]^2]
Print[bs]
f0[x_]=If[x+a[8]<-2.33810, 0, n0 AiryAi[x+a[8]]  ]
For[i=1,i<16,i++,f[i,x_]=n[i]*AiryAi[x+a[i]] ]
p[t_,x_]=Sum[b[i]*f[i,x]*Exp[I(a[i]-a[8])t],{i,1,15}]
Plot[{Abs[f0[x]]^2,Abs[p[0,x]]^2},{x,80,130},PlotRange->{0,.65},PlotStyle->{{RGBColor[0,0,0]},{RGBColor[1,0,0]}} ]
Plot[Abs[p[0,x]]^2,{x,80,130},PlotRange->{0,.65}]
Plot[Abs[p[1,x]]^2,{x,80,130},PlotRange->{0,.65}]
Plot[Abs[p[2,x]]^2,{x,80,130},PlotRange->{0,.65}]
Plot[Abs[p[3,x]]^2,{x,80,130},PlotRange->{0,.65}]
Plot[Abs[p[4,x]]^2,{x,80,130},PlotRange->{0,.65}]
temp1[x_]=Sum[Abs[p[4,x+i]]^2,{i,0,130,5}]
norm=NIntegrate[temp1[x],{x,0,5}]
Sum[Abs[b[i]]^2,{i,1,15}]
temp2[t_,x_]=Sum[(x+i)*Abs[p[t,x+i]]^2,{i,0,130,5}]
avgx[t_]=NIntegrate[temp2[t,x],{x,0,5}]/norm
Plot[avgx[t],{t,0,8}]
----
NIntegrate[p[0,x]^2,{x,0,130}]
(0.954693 vs sum of squares of bs=0.954779)
NIntegrate[p[0,x]^2 x,{x,0,130}]/%
(121.253)
NIntegrate[p[0,x]^2 (x-121.253)^2,{x,0,130}]/%%
(14.8829, Delta x=3.8578)
junk[x_]=D[p[0,x],x]
NIntegrate[junk[x]^2,{x,0,130}]/0.954693
(0.815433=<p^2>, Delta p =.90301, Delta v=2 Delta p = 1.8060)
--
NIntegrate[p[0,x]^2 ,{x,116,126}]
(0.945172)
NIntegrate[p[0,x]^2 x,{x,116,126}]/%
(121.499)
NIntegrate[p[0,x]^2 (x-121.5)^2,{x,116,126}]/%%
(0.762371=<x^2>, Delta x= 0.873138)
NIntegrate[junk[x]^2,{x,116,126}]/0.945172
( 0.56104=<p^2>, Delta p=0.749026)
--
NIntegrate[f0[x]^2,{x,116,126}]
(0.999999)
NIntegrate[f0[x]^2 x,{x,116,126}]/%
(121.446)
NIntegrate[f0[x]^2 (x-121.446)^2,{x,116,126}]/%%
(0.485908=<x^2>, Delta x=0.697071)
NIntegrate[f0'[x]^2,{x,116,126}]/0.999999
(0.783734=<p^2>, Delta p=0.885288)

---
<<Statistics`ContinuousDistributions`
s[t_]=Sqrt[(.873141)^2+(1.5 t)^2]
m[t_]=121.5-t^2
pc[t_,x_]=PDF[NormalDistribution[m[t],s[t]],x]
Plot[{pc[0,x],pc[1,x],pc[2,x],pc[4,x]},{x,80,130},PlotRange->{0,.8}]]

-----
f[x_,a_]=x(a^2-x^2)
Integrate[f[x,a] (-D[D[f[x,a],x],x]+x f[x,a]),{x,0,a}]
Integrate[D[f[x,a],x]^2+x f[x,a]^2,{x,0,a}]
(both: (4/5)a^5+a^8/24 )
Integrate[f[x,a]^2,{x,0,a}]
(8a^7/105)
Simplify[%%/%]
[(21/2)a^{-2}+(35/64)a)
fn[x_,a_]=If[x>0 && x<a ,x(a^2-x^2)/Sqrt[(8a^7/105)],0]
a/: Positive[a]=True
Integrate[g[x,a] (-D[D[g[x,a],x],x]+x g[x,a]),{x,0,Infinity}]
Integrate[g[x,a]^2,{x,0,Infinity}]
Simplify[%%/%]
[(3/2)a^(-1)+a^2)
 % /. a->(3/4)^(1/3)
2.47645
gn[x_,a_]=x Exp[-a x] Sqrt[(4 a^3)]
n0=1/Sqrt[NIntegrate[AiryAi[x]^2,{x,-2.33810,6}]]
Plot[{fn[x,4(3/5)^(1/3)],gn[x,(3/4)^(1/3)],n0 AiryAi[x-2.3381]},
{x,0,6},PlotStyle->{{RGBColor[0,0,1]},{RGBColor[1,0,0]},{RGBColor[0,0,0]}} ]
----
FindRoot[Sqrt[10+x]-AiryAi'[x]/AiryAi[x], {x,-2.3}]
For[i=1,i<8,i++,b[i]=x /. FindRoot[Sqrt[10+x]-AiryAi'[x]/AiryAi[x], {x,a[i]+.5}]]
----
For[i=1,i<8,i++,c[i]=x /. FindRoot[AiryAi'[x]/AiryAi[x]+3 AiryAi'[x/9]/AiryAi[x/9], {x,a[i]+.5}]]
----
For[i=1,i<17,i=2 i,xavg[i]=NIntegrate[f[i,x]^2 x,{x,0,-a[i]+5}]]
For[i=1,i<17,i=2 i,x2[i]=NIntegrate[f[i,x]^2 (x-xavg[i])^2,{x,0,-a[i]+5}]]
For[i=1,i<17,i=2 i,fp[i,x_]=D[f[i,x],x]]
For[i=1,i<17,i=2 i,p2[i]=NIntegrate[fp[i,x]^2 ,{x,0,-a[i]+5}]]
For[i=1,i<17,i=2 i,Print[{i,Sqrt[x2[i]] Sqrt[p2[i]] } ] ]
---

FindRoot[Sqrt[x-a1]Cot[a1 Sqrt[x-a1]] -AiryAi'[a1-x]/AiryAi[a1-x], {x,2.3}]
For[i=1,i<8,i++,ea1[i]=x /. FindRoot[Sqrt[x-a1]Cot[a1 Sqrt[x-a1]] -AiryAi'[a1-x]/AiryAi[a1-x], {x,-a[i]+.3}]]
For[i=1,i<8,i++,ea2[i]=-a[i]+NIntegrate[(a1-x)f[i,x]^2,{x,0,a1}]]
ex1[0]=-a[1]
For[i=1,i<17,i++,a1=.3 i; ex1[i]=x /. FindRoot[Sqrt[x-a1]Cot[a1 Sqrt[x-a1]] -AiryAi'[a1-x]/AiryAi[a1-x], {x,ex1[i-1]}]]
ePT1[0]=-a[1]
For[i=1,i<17,i++,a1=.3 i; ePT1[i]=-a[1]+NIntegrate[(a1-x)f[1,x]^2,{x,0,a1}]]
ex2[0]=-a[2]
For[i=1,i<17,i++,a1=.3 i; ex2[i]=x /. FindRoot[Sqrt[x-a1]Cot[a1 Sqrt[x-a1]] -AiryAi'[a1-x]/AiryAi[a1-x], {x,ex2[i-1]+.1}]]
ePT1[0]=-a[2]
For[i=1,i<17,i++,a1=.3 i; ePT2[i]=-a[2]+NIntegrate[(a1-x)f[2,x]^2,{x,0,a1}]]
tex1=Table[{.3 i,ex1[i]},{i,0,16}]
tex2=Table[{.3 i,ex2[i]},{i,0,16}]
tPT1=Table[{.3 i,ePT1[i]},{i,0,16}]
tPT2=Table[{.3 i,ePT2[i]},{i,0,16}]
Show[ListPlot[tPT1,PlotJoined->True,PlotStyle->{RGBColor[1,0,0]},PlotRange->{0,6.5}],
ListPlot[tPT2,PlotJoined->True,PlotStyle->{RGBColor[1,0,0]},PlotRange->{0,6.5}],
ListPlot[tex1,PlotJoined->True,PlotRange->{0,6.5}],
ListPlot[tex2,PlotJoined->True,PlotRange->{0,6.5}] ]
