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=, 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=, Delta x= 0.873138) NIntegrate[junk[x]^2,{x,116,126}]/0.945172 ( 0.56104=, 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=, Delta x=0.697071) NIntegrate[f0'[x]^2,{x,116,126}]/0.999999 (0.783734=, Delta p=0.885288) --- <{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(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}] ]