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 quadrupole oscillation mode of an axially symmetric BEC. It uses the auxiliary file "BEC-Frequency-Shift-Axial-Q-Init.m", produced by the notebook "BEC-Frequency-Shift-Axial-Q-Init.nb".

In order  to study a nonlinear BEC dynamics analytically, we start from the variationally-obtained system of nonlinear differential equations
   D[u[t], {t, 2}] + u[t] - 1/u[t]^3 - p/(u[t]^3 v[t]) - q Cos[Ω t]/(u[t]^3 v[t]) == 0,
   D[v[t], {t, 2}] + λ^2 v[t] - 1/v[t]^3 - p/(u[t]^2 v[t]^2) - q Cos[Ω t]/(u[t]^2 v[t]^2) == 0.

We next choose the dimensionless parameters that desribe the condensate - interaction strength p and trap anisotropy λ:

p=1;
λ=3/10;

From the linear stability analysis, we obtain equilibrium sizes u0 (in the radial) and v0 (in the z direction), and also frequencies of the quadrupole and breathing mode (ωQ0 and ωB0, respectively):

Clear[u0,v0];
eqsolution=NSolve[{-(1/u0^3)+u0-p/(u0^3 v0)==0,-(1/v0^3)-p/(u0^2 v0^2)+v0 λ^2==0},{u0,v0},WorkingPrecision->50];
l=Length[eqsolution];
For[i=1,i≤l,i++,{u0,v0}={u0,v0}/.eqsolution[[i]];
If[(u0∈Reals)&&(v0∈Reals),If[(u0>0)&&(v0>0),Print[u0," ",v0];mem=i;]];
Clear[u0,v0];];
{u0,v0}={u0,v0}/.eqsolution[[mem]];
p23=p/(4 u0^2 v0^3);
p32=p/(4 u0^3 v0^2);
ωQ0=Sqrt[2] Sqrt[(1+λ^2-p23)-Sqrt[(1-λ^2+p23)^2+8 p32^2]];
ωB0=Sqrt[2] Sqrt[(1+λ^2-p23)+Sqrt[(1-λ^2+p23)^2+8 p32^2]];
Print[N[p]," ",N[u0]," ",N[v0]," ",N[ωQ0]," ",N[ωB0]];

In the next step we find numerical constants necessery for calculation of the shift of the quadrupole mode:

m11=4;
m12= p/(u0^3 v0^2);
m21=2m12;
m22=λ^2+3/v0^4+2p/(u0^2  v0^3);
EE=Sqrt[(m22-m11)^2+4m12 m21];
c1=((m22-m11)-EE)/(2m12);
c2=((m22-m11)+EE)/(2 m12);
Ω1= ωQ0;
Ω2= ωB0;

We import definitions of constants Aq and Cq from the auxiliary file "BEC-Frequency-Shift-Axial-Q-Init.m":

<<"BEC-Frequency-Shift-Axial-Q-Init.m";

aux=0;
For[i=1,i≤13,i++,aux+=Cq[[i]]];

Clear[Cq];Cq=aux;

Finally, the frequency shift is calculated as:

QMFrequencyShift[Ω_,q_]=Simplify[q^2Expand[(Cq/(2 Ω1 Aq))]]

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

Plot[{ωQ0,ωQ0+QMFrequencyShift[Ω,2/10], ωQ0+QMFrequencyShift[Ω,8/10]},{Ω,0,3},PlotRange-> {ωQ0-ωQ0/10,ωQ0+ωQ0/10},PlotStyle-> {{Dashed,Thick,Green},{Thick,Blue},{Dashed,Thick,Red}}]