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.

This initialization part of the code is used for the calculation of the coefficients Ab and Cb, which are used as inputs for obtaining the frequency shift of the breathing oscillation mode of an axially symmetric BEC. It produces the auxiliary file "BEC-Frequency-Shift-Axial-B-Init.m", where the coefficients are stored for further use by the main part of the code, notebook "BEC-Frequency-Shift-Axial-B.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 use the following perturbative expansions in q in order to obtain a system of linear differential equations:

u[t]=u0+q u1[t]+q^2 u2[t]+q^3 u3[t];
v[t]=v0+q v1[t]+q^2 v2[t]+q^3 v3[t];
Equation1= 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]);
Equation2= 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);

From the zeroth order, we obtain the equilibrium condensate size u0 (in the radial direction) and v0 (in the z direction):

Normal[Series[Equation1,{q,0,0}]]

Normal[Series[Equation2,{q,0,0}]]

First order expansion yields a system of equations for u1[t] and v1[t]:

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

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

It turns out, that in each order n > 0, we have the system of equations of the general form 
un''[s] + m11 un[s] + m12 vn[s] + fnu[s] = 0,
vn''[s] + m21 un[s] + m22 vn[s] + fnv[s] = 0.

f1u[t]=-((Cos[t Ω])/((u0^3) v0));
f1v[t]=-((Cos[t Ω])/((u0^2) (v0^2)));

In the second order of the perturbation theory we obtain the following equations for u2[t] and v2[t]:

Expand[Coefficient[Normal[Series[Equation1,{q,0,2}]],q^2]]

Expand[Coefficient[Normal[Series[Equation2,{q,0,2}]],q^2]]

f2u[t]=((3 Cos[t Ω] u1[t])/((u0^4) v0))-((6 ((u1[t])^2))/(u0^5))-((6 p ((u1[t])^2))/((u0^5) v0))+((Cos[t Ω] v1[t])/((u0^3) (v0^2)))-((3 p u1[t] v1[t])/((u0^4) (v0^2)))-((p ((v1[t])^2))/((u0^3) (v0^3)));
f2v[t]=((2 Cos[t Ω] u1[t])/((u0^3) (v0^2)))-((3 p ((u1[t])^2))/((u0^4) (v0^2)))+((2 Cos[t Ω] v1[t])/((u0^2) (v0^3)))-((4 p u1[t] v1[t])/((u0^3) (v0^3)))-((6 ((v1[t])^2))/(v0^5))-((3 p ((v1[t])^2))/((u0^2) (v0^4)));

In the third order of the perturbation theory we find the following equations for u3[t] and v3[t], which contain f3u[t] and f3v[t]:

f3u=Expand[Coefficient[Normal[Series[Equation1,{q,0,3}]],q^3]]

f3upart[1,t]=-((6 Cos[t Ω] ((u1[t])^2))/((u0^5) v0));
f3upart[2,t]=((10 ((u1[t])^3))/(u0^6))+((10 p ((u1[t])^3))/((u0^6) v0));
f3upart[3,t]=-((3 Cos[t Ω] u1[t] v1[t])/((u0^4) (v0^2)));
f3upart[4,t]=((6 p ((u1[t])^2) v1[t])/((u0^5) (v0^2)));
f3upart[5,t]=-((Cos[t Ω] ((v1[t])^2))/((u0^3) (v0^3)));
f3upart[6,t]=((3 p u1[t] ((v1[t])^2))/((u0^4) (v0^3)));
f3upart[7,t]=((p ((v1[t])^3))/((u0^3) (v0^4)));
(* u2 terms *)
f3upart[8,t]= ((3 Cos[t Ω] )/((u0^4) v0));
f3upart[9,t]= -((12 u1[t] )/(u0^5))-((12 p u1[t] )/((u0^5) v0));
f3upart[10,t]= -((3 p  v1[t])/((u0^4) (v0^2)));
(* v2 terms *)
f3upart[11,t]= ((Cos[t Ω] )/((u0^3) (v0^2)));
f3upart[12,t]= -((3 p u1[t])/((u0^4) (v0^2)));
f3upart[13,t]= -((2 p v1[t] )/((u0^3) (v0^3)));

f3v=Expand[Coefficient[Normal[Series[Equation2,{q,0,3}]],q^3]]

f3vpart[1,t]=-((3 Cos[t Ω] ((u1[t])^2))/((u0^4) (v0^2)));
f3vpart[2,t]=((4 p ((u1[t])^3))/((u0^5) (v0^2)));
f3vpart[3,t]=-((4 Cos[t Ω] u1[t] v1[t])/((u0^3) (v0^3)));
f3vpart[4,t]=((6 p ((u1[t])^2) v1[t])/((u0^4) (v0^3)));
f3vpart[5,t]=-((3 Cos[t Ω] ((v1[t])^2))/((u0^2) (v0^4)));
f3vpart[6,t]=((6 p u1[t] ((v1[t])^2))/((u0^3) (v0^4)));
f3vpart[7,t]=((10 ((v1[t])^3))/(v0^6))+((4 p ((v1[t])^3))/((u0^2) (v0^5)));
(* u2 terms *)
f3vpart[8,t]= ((2 Cos[t Ω] )/((u0^3) (v0^2)));
f3vpart[9,t]=-((6 p u1[t] )/((u0^4) (v0^2)));
f3vpart[10,t]=-((4 p  v1[t])/((u0^3) (v0^3)));
(* v2 terms *)
f3vpart[11,t]=((2 Cos[t Ω] )/((u0^2) (v0^3)));
f3vpart[12,t]=-((4 p u1[t] )/((u0^3) (v0^3)));
f3vpart[13,t]=-((12 v1[t] )/(v0^5))-((6 p v1[t] )/((u0^2) (v0^4)));

In order to solve previous equations, we start from the first order and, instead of the coupled equations for u1[t] and v1[t], use the decoupled equations for x1[t] and y1[t]:

f1x[t]=( c2 f1u[t]-f1v[t])/(c2-c1);
f1y[t]=-( c1 f1u[t]-f1v[t])/(c2-c1);

We solve explicitely for x1[t] and y1[t]:

x1[t_]=x1[t]/.Simplify[Expand[TrigReduce[DSolve[{f1x[t]+Ω1^2x1[t]+x1''[t]==0,x1[0]==0,x1'[0]==0},x1[t],t]]]][[1]]

y1[t]=y1[t]/.Simplify[Expand[TrigReduce[DSolve[{f1y[t]+Ω2^2y1[t]+ y1''[t]==0,y1[0]==0,y1'[0]==0},y1[t],t]]]][[1]]

Now, we have u1[t] and v1[t]:

u1[t]=x1[t]+y1[t];
v1[t]=c1 x1[t]+c2 y1[t];

In the second order, we transform the coupled equations for u2[t] and v2[t] into the decoupled equations for x2[t] and y2[t]:

f2x[t]=( c2 f2u[t]-f2v[t])/(c2-c1);
f2y[t]=-( c1 f2u[t]-f2v[t])/(c2-c1);

For practical reasons, we represent f2x[t] and f2y[t] as scalar products f2xcoeff × functions and f2ycoeff × functions, where:

f2xcoeff={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
f2ycoeff={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
functions={Cos[t Ω1+t Ω2],Sin[t Ω1+t Ω2],Cos[t Ω+t Ω2],Sin[t Ω+t Ω2],Cos[t Ω1-t Ω2],Sin[t Ω1-t Ω2],Cos[t Ω-t Ω2],Sin[t Ω-t Ω2],Cos[2t Ω2], Sin[2t Ω2],Cos[2 t Ω1], Sin[2t Ω1],Cos[2 t Ω], Sin[2t Ω],Cos[t Ω+t Ω1],Sin[t Ω+t Ω1],Cos[t Ω-t Ω1],Sin[t Ω-t Ω1]};

Length[functions]

f2xx=f2x[t];
f2yy=f2y[t];
imax=Length[functions];
For[i=1,i≤imax,i++,
(*Print[i];*)
f2xx=TrigReduce[f2xx];
f2yy=TrigReduce[f2yy];
f2xcoeff[[i]]=Coefficient[f2xx,functions[[i]]];
f2ycoeff[[i]]=Coefficient[f2yy,functions[[i]]];
f2xx=Expand[f2xx-f2xcoeff[[i]] functions[[i]]];
f2yy=Expand[f2yy-f2ycoeff[[i]]functions[[i]]];
];
constf2y=f2yy;
constf2x=f2xx;

Now, we solve explicitely for x2[t] and y2[t]:

functions2={1,Cos[Ω1 t], Cos[Ω2 t],Cos[t Ω1+t Ω2],Cos[t Ω+t Ω2],Cos[t Ω1-t Ω2],Cos[t Ω-t Ω2],Cos[2 t Ω2],Cos[2 t Ω1],Cos[2 t Ω],Cos[t Ω+t Ω1],Cos[t Ω-t Ω1]};
frequencies={0,Ω1, Ω2, Ω1+Ω2, Ω+Ω2,Ω1-Ω2,Ω-Ω2,2Ω2,2Ω1,2Ω,Ω+Ω1,Ω-Ω1};
coefficientsx2={0,0,0,0,0,0,0,0,0,0,0,0};
coefficientsy2={0,0,0,0,0,0,0,0,0,0,0,0};
coefficientsx2[[1]]=-constf2x/Ω1^2;
coefficientsy2[[1]]=-constf2y/Ω2^2;
coefficientsx2[[2]]-=coefficientsx2[[1]];
coefficientsy2[[3]]-=coefficientsy2[[1]];
For[i=4,i≤12,i++,
coefficientsx2[[i]]=-(1/(Ω1^2-frequencies[[i]]^2)) f2xcoeff[[Position[functions,functions2[[i]]][[1]][[1]]]];
coefficientsy2[[i]]=-(1/(Ω2^2-frequencies[[i]]^2)) f2ycoeff[[Position[functions,functions2[[i]]][[1]][[1]]]];
coefficientsx2[[2]]-=coefficientsx2[[i]];
coefficientsy2[[3]]-=coefficientsy2[[i]];
];

Finally, we solve explicitely for u2[t] and v2[t]:

coefficientsu2={0,0,0,0,0,0,0,0,0,0,0,0};
coefficientsv2={0,0,0,0,0,0,0,0,0,0,0,0};
For[i=1,i≤12,i++,
coefficientsu2[[i]]=Expand[Simplify[coefficientsx2[[i]]+coefficientsy2[[i]]]];
coefficientsv2[[i]]=Expand[Simplify[c1 coefficientsx2[[i]]+c2 coefficientsy2[[i]]]];
];

A term Cos[Ω2 t]] appears in the equation for x3[t]. We calculate the coefficient Cb  in front of it, and for practical reasons keep it in the split form:

Cb={0,0,0,0,0,0,0,0,0,0,0,0,0};
For[j=1,j≤ 7,j++,
Cb[[j]]+=Coefficient[TrigReduce[c1 f3upart[j,t]-f3vpart[j,t]],Cos[Ω2 t]]/(c1-c2);
Print[j];
];

For[j=8,j≤ 10,j++,
Print[j];
For[k=1,k≤ 12,k++,
Cb[[j]]+=coefficientsu2[[k]](c1 Coefficient[TrigReduce[f3upart[j,t] functions2[[k]]],Cos[Ω2 t]]-
Coefficient[TrigReduce[f3vpart[j,t] functions2[[k]]],Cos[Ω2 t]])/(c1-c2);
];
];

For[j=11,j≤ 13,j++,
Print[j];
For[k=1,k≤ 12,k++,
Cb[[j]]+=coefficientsv2[[k]](c1 Coefficient[TrigReduce[f3upart[j,t] functions2[[k]]],Cos[Ω2 t]]-
Coefficient[TrigReduce[f3vpart[j,t] functions2[[k]]],Cos[Ω2 t]])/(c1-c2);
];
];

For further calculations, we save the results: the coefficient in front of the Cos[Ω2 t] term in the first order solution y1[t], i.e. Ab, as well as the coefficient in front of Cos[Ω2 t] in the equation for y3[t], i.e. Cb:

Ab=Coefficient[y1[t],Cos[t Ω2]]; 
Definition[Ab]>>>"BEC-Frequency-Shift-Axial-B-Init.m";
Cb=Expand[Cb];
Definition[Cb]>>>"BEC-Frequency-Shift-Axial-B-Init.m";