This code is written by Ivana Vidanovic (ivanavi@ipb.ac.rs) from the Scientific Computing Laboratory of the Institute of Physics Belgrade, http://www.scl.rs/  for the paper "Nonlinear BEC Dynamics Induced by a Harmonic Modulation of the s-wave Scaterring Length" by I. Vidanovic, A. Balaz, H. Al-Jibbouri, and A. Pelster.

The code is used for calculation of the frequency shift of the collective oscillation mode of a spherically symmetric BEC.

In order  to study a nonlinear BEC dynamics analytically, we start from the variationally-obtained nonlinear differential equation
   D[u[t], {t, 2}] + u[t] - 1/u[t]^3 - p/u[t]^4 - q/u[t]^4 Cos[Ω  t] == 0.
We apply the Poincare-Lindstedt method and first preform a rescaling of time s=ω t:

Equation=ω^2 D[u[s],{s,2}]+u[s]-1/u[s]^3-p/u[s]^4-q/u[s]^4 Cos[Ω  s/ωω];

(* Explanation: ω=ωω, but we do not want to expand Cos[Ω s/ω'] with respect to q. *)

We use the following perturbative expansions in q in order to obtain a system of linear differential equations:

u[s_]=u0+q u1[s]+q^2u2[s]+q^3 u3[s];
ω = ω0+q ω1+q^2 ω2+q^3 ω3;

From the zeroth order, we obtain equilibrium condensate size u0:

Coefficient[Expand[Normal[Series[Equation,{q,0,3}]]],q,0]

First order expansion yields:

Coefficient[Expand[Normal[Series[Equation,{q,0,3}]]],q,1]

In the second order, we have:

Coefficient[Expand[Normal[Series[Equation,{q,0,3}]]],q,2]

and in the third order, we obtain:

Coefficient[Expand[Normal[Series[Equation,{q,0,3}]]],q,3]

In each order we have an equation of the general form ω0^2 un''[s]+ω0^2 un[s]+fn[s]=0, n>0, where:

f1[s]=-((Cos[((s Ω)/(ωω))])/(u0^4));
f2[s]=((4 Cos[((s Ω)/(ωω))] u1[s])/(u0^5))-((10 p ((u1[s])^2))/(u0^6))-((6 ((u1[s])^2))/(u0^5))+2 ω0 ω1 u1''[s];
f3[s]=-((10 Cos[((s Ω)/(ωω))] ((u1[s])^2))/(u0^6))+((20 p ((u1[s])^3))/(u0^7))+((10 ((u1[s])^3))/(u0^6))+((4 Cos[((s Ω)/(ωω))] u2[s])/(u0^5))-((20 p u1[s] u2[s])/(u0^6))-((12 u1[s] u2[s])/(u0^5))+((ω1)^2) u1''[s]+2 ω0 ω2 u1''[s]+2 ω0 ω1 u2''[s];

We solve for u1[s]:

Clear[u1];u1[s_]=u1[s]/.Simplify[Expand[TrigReduce[DSolve[{f1[s]+ω0^2u1[s]+((ω0)^2) u1''[s]==0,u1[0]==0,u1'[0]==0},u1[s],s]]]][[1]]

and for u2[s]:

Clear[u2];u2[s_]=u2[s]/.TrigReduce[DSolve[{f2[s]+ω0^2u2[s]+ω0^2 u2''[s]==0,u2[0]==0,u2'[0]==0},u2[s],s]][[1]]

From the expression for u2[s], we conclude that:

ω1=0;

which leads to the second-order result without secular terms:

u2[s_]=TrigReduce[Simplify[u2[s]]]

In the third order of the perturbation theory, a term Cos[s] appears in f3[s] with the following coefficient:

Coeff3=Simplify[Coefficient[TrigReduce[Expand[f3[s]]],Cos[s]]]

In order to elliminate secular terms from the solution, the value of ω2 is determined so that the coefficient in front of Cos[s] vanishes:

FrequencyShift[Ω_,q_]=q^2Simplify[ω2/.Solve[Coeff3==0,ω2][[1]]/.{ωω-> ω0}]

An example of the frequency shift calculation for p = 1:

p=1;
u0=u/.NSolve[u-1/u^3-p/u^4==0,u][[1]];
ω0=Sqrt[1+3/u0^4+4 p/u0^5];

Plot of the breathing mode frequency versus the external modulation frequency Ω for q = 0.4 (blue) and q = 0.8 (red):

Plot[{ω0,ω0 +FrequencyShift[Ω,4/10],ω0 +FrequencyShift[Ω,8/10]},{Ω,0.9,8},PlotRange->{ω0-ω0/10,ω0+ω0/10},PlotStyle->{{Black,Thick},{Blue,Thick,Dashed},{Red,Thick,Dashed}}]