Free fall

Dynamics

Exact solution

Clear[eigen];
\(\(eigen[m_Integer] := \(eigen[m] = 
          Block[{x0, xy}, 
            x0\  = \ 
              N[\((\((2  m - 1\/2)\) \(π\ 3\)\/4)\)\^\(2/3\)]; 
xy /. 
              FindRoot[
                AiryAi[xy] == 
                  0, {xy, \(-x0\)}]
						]\);\)\)

Clear[ϕ];
ϕ[m_Integer, x_]  :=  AiryAi[x + eigen[m]]

Clear[z];
z[m_Integer] :=  z[m]  =  
          Sqrt[ NIntegrate[ 
              Abs[ϕ[m, x] ]^2, {x, 0, Infinity}]];

SetOptions[NIntegrate, AccuracyGoal -> 10, 
        MaxRecursion -> 9];
Clear[newpsi];
newpsi[x0_:  7, σ_:  1, 
            mmax_:  15]  /;  (NumberQ[x0] && NumberQ[σ] && 
              IntegerQ[
                15]) := (
Clear[ψ]; 
ψ[x_] :=  
            Exp[-(x - 
                        x0)^2/(2  σ^2)]; 
norm  =  
            Sqrt[ NIntegrate[
                Abs[ψ[x] ]^2, {x, 0, 
                  Infinity}]]; 
ψ[x]/
            norm; 
 (**) 
Clear[
            c]; 
c[m_Integer] := c[m]  =  
              NIntegrate[ 
                Conjugate[ϕ[m, x]/z[m]] ψ[x]/norm, {x, 0, 
                  Infinity}]; 
clist  =  
            Table[c[m], {m, 1, mmax}]; 
taulist  =  
            Table[2  Pi/eigen[m], {m, 1, 
                mmax}]; 
 (* MultipleListPlot[clist, 
                taulist]; *) 
 (**) 
ψ[x_, t_]  :=  
            Sum [ c[m] Exp[-I t eigen[m] ] ϕ[m, x] /
                  z[m] , {m, 1, 
                mmax}]; 
ψmax  =  ψ[x0, 
              0]; 
 (**) 
Clear[
            err]; 
err[n_]  :=  err[n]  =  
              err[n - 1]  -  c[n]^2; 
err[0] = 
            1; 
err[
            mmax]
);

    
newpsi[]

Clear[xmean, xstdev2, xstdev];
xmean[t_?NumberQ]  :=  xmean[t] = 
          NIntegrate[ 
            x   Abs[ψ[x, t]]^2, {x, 0, 
              Infinity}];
xstdev2[t_?NumberQ]  :=  
        NIntegrate[ x^2   Abs[ψ[x, t]]^2, {x, 0, Infinity}] -  
          xmean[t]^2;
xstdev[t_] :=  xstdev[t] = Sqrt[xstdev2[t]];

Clear[plot, tplot];
plot[t_, xmax_:  15, nx_:  50, opts___]  :=  
        Plot[Abs[ψ[x, t]]^2, {x, 0, xmax}, PlotPoints -> nx, 
          PlotRange -> {{0, xmax}, {0,  ψmax}} , 
          opts];

    
dplot[t_, xmax_:  15, nx_:  50, 
          opts___]  :=  (xmeant = xmean[t]; xdev = xstdev[t]; 
          Show[plot[t, xmax, nx, DisplayFunction -> Identity, 
              PlotStyle -> pstyle], 
Graphics[{lstyle, 
                Line[{{xmeant, 0}, {xmeant, ψmax}}], 
                Line[{{xmeant - xdev, ψmax}, {xmeant + 
                        xdev, ψmax}}]}], 
            opts]);

    
pstyle = RGBColor[0. , 0. , 0. ]; 
    lstyle = RGBColor[1. , 0, 
        0];

tplot[tmax_:  40, dt_:  .5, xmax_:  15, 
        nx_:  50] :=  
Table[
        plot[t, xmax, nx, DisplayFunction -> Identity], {t, 0, tmax, 
          dt}]

dplot[10. , 15, 50, 
        DisplayFunction -> $DisplayFunction];

plotall = tplot[];

pstyle = GrayScale[1]; lstyle = Hue[1]; 
    Show[plotall, DisplayFunction -> $DisplayFunction]

Clear[dfun];
dfun[n_] := $DisplayFunction;
Do[
        Show[plotall[[ii]], DisplayFunction -> dfun[ii]], {ii, 
          Length[plotall]}];
  
SelectionMove[EvaluationNotebook[], After, 
        EvaluationCell];
SelectionMove[EvaluationNotebook[], Next, CellGroup];
FrontEndTokenExecute["OpenCloseGroup"];
\(	\)