Clear["Global`*"]

(* Number of generators, i.e. the dimension of the algebra: 
for su(2) Lie algebra we set basisNo = 3
for su(3) Lie algebra we set basisNo = 8
*)
basisNo = 3;

(* Highest Fourier component to be taken into account:
   all Fourier components with indices m such that |m| > fourierCompNo remain zero throughout the flow
*)
fourierCompNo = 2;

(* Power series expansion order:
  we obtain the inverse-frequency expansion of the effective Hamiltonians
with terms up to, and including, those proportional to omega^{-orderNo}
*)
orderNo =5;

(* Expansion of the micromotion operator:
by default, the micromotion operator is obtained to the same order as the effective Hamiltonian, thus we magnusTermNo = orderNo
win order to limit the expansion of the micromotion to fewer terms, set magnusTermNo to a lower value
*)
magnusTermNo =orderNo;

(* Specify which parameters are complex 
	by default none, i.e. not needed in our examples
*)
complexParameters = {};

(* Specify initial conditions at the begining of the flow, that is for the flow variable s = 0
Only non-zero-valued Fourier components need to be listed as the rest will be included automatically below.
Use the format: x_j [n] ==  some function of parameters and possibly of time t
Here j is the index of the element of the algebra, that is x_j is the coefficient multiplying the generator G_j
In the brackets, n is the index of the Fourier component -- include only zero and positive indices, as the rest will be generated automatically.
 For example, the model of Eqs. (60) in the manuscript is described 
initCond = {(Subscript[x, 1])[0]== g[t]*Cos[ϕ],
	(Subscript[x, 2])[0]== g[t]*Sin[ϕ],
	(Subscript[x, 3])[0]== Δ/2,
	(Subscript[x, 1])[2]== g[t]/2*Cos[ϕ]+I*g[t]/2*Sin[ϕ],
	(Subscript[x, 2])[2]== -g[t]/2*Sin[ϕ]+I*g[t]/2*Cos[ϕ]} ;
Please note that writing g[t] implies that the parameter g is time dependent, and the resulting expansions will promptly include its time derivatives.
In case of constant coupling, please replace g[t] with just g. Likewise, in the above example phi is a constant. 
*)
initCond = {(Subscript[x, 1])[0]== g[t]*Cos[ϕ],
		(Subscript[x, 2])[0]== g[t]*Sin[ϕ],
		(Subscript[x, 3])[0]== Δ/2,
		(Subscript[x, 1])[2]== g[t]/2*Cos[ϕ]+I*g[t]/2*Sin[ϕ],
		(Subscript[x, 2])[2]== -g[t]/2*Sin[ϕ]+I*g[t]/2*Cos[ϕ]};

(* Commutator -- implementation of the Lie algebra
	The function receives two lists of coefficients, p_i and q_i, that multiply the generators G_i
	and returns the commutator [p, q] as the corresponding list of coefficients c: \sum_k c_k G_k = [p_i G_i, q_j G_j]
The example below implements the usual su(2) algebra. 
*)
Comm[p_,q_]:=
{2*I*p[[2]]*q[[3]]-2*I*q[[2]]*p[[3]],
2*I*p[[3]]*q[[1]]-2*I*q[[3]]*p[[1]],
2*I*p[[1]]*q[[2]]-2*I*q[[1]]*p[[2]]}

(* Fill the table of initial conditions by including the specified ones and generate zero-valued initial conditions. *)
initCondTemp ={};
Do[AppendTo[initCondTemp,i[[1]][[0]][i[[1]][[1]],0,0,t] == i[[2]]],{i,initCond}]
initCondAll = Flatten[Table[(Subscript[x, k])[i,j,0,t]==0,{i,0,fourierCompNo},{j,0,orderNo},{k,1,basisNo}]];
Do[posRep =Position[initCondAll[[;;,1]],initCondTemp[[i,1]]];
initCondAll[[posRep[[1]]]]=initCondTemp[[i]];
,{i, Dimensions[initCondTemp][[1]]}]

(* Abstract functions to hold the solutions. *)
solFunctionsPower = Flatten[Table[(Subscript[x, k])[i,j,s,t],{i,0,fourierCompNo},{j,0,orderNo},{k,1,basisNo}]];
solFunctions = Flatten[Table[{(Subscript[c, k])[i,s]},{i,0,fourierCompNo},{k,1,basisNo}]];

(* Power expansion of terms. *)
xx[k_][m_,s]:= Sum[ϵ^i*(Subscript[x, k])[m,i,s,t],{i,0,orderNo}];

(* General operator form. *)
H[m_,s] := Table[xx[i][m,s],{i,1,basisNo}]/; m≥0;
H[m_,s] := Conjugate[Table[xx[i][-m,s] ,{i,1,basisNo}]] /; m<0;

(* Flow equation for the zeroth Fourier component, Eq. (54) of the manuscript. *)
flowEqs = Thread[D[H[0,s],s] == 2*ϵ*Sum[Comm[H[m,s],H[-m,s]],{m,1,fourierCompNo}]];

(* Flow equations for other Fourier components, Eq. (55) of the manuscript. *)
Do[
eqTemp = 
 Thread[D[H[n,s],s] ==-n* H[n,s]+ϵ*Comm[H[n,s],H[0,s]]+I*ϵ*D[H[n,s],t]+2*ϵ*Sum[Comm[H[n+l,s],H[-l,s]],{l,1,fourierCompNo-n}]]; 
flowEqs=Join[flowEqs,eqTemp];
,{n,1,fourierCompNo}
]

(* Begin solving equations. *)
(* Solve for zeroth term. *)
initConditions  = Pick[initCondAll,initCondAll[[;;,1]][[;;,2]],0];
resFunctions = Pick[solFunctionsPower ,solFunctionsPower [[;;,2]],0];
sol = DSolve[Join[D[flowEqs,{ϵ,0}]/. ϵ -> 0,initConditions],resFunctions,s];
updatedSol =  Join[sol[[1]],D[sol[[1]],t]];
finalRes = Pick[solFunctionsPower ,solFunctionsPower [[;;,2]],0] /. sol[[1]];

(*Solve rest of terms*)
Do [
Print["calculating order:",i];
initConditions  = Pick[initCondAll,initCondAll[[;;,1]][[;;,2]],i];
resFunctions = Pick[solFunctionsPower ,solFunctionsPower [[;;,2]],i];
flowEqsUpdated :=D[ComplexExpand[flowEqs/. updatedSol,complexParameters] ,{ϵ,i}]/. ϵ -> 0;
sol = DSolve[Join[flowEqsUpdated ,initConditions],resFunctions,s];
updatedSol =  Join[updatedSol ,sol[[1]],D[sol[[1]],t]];
finalRes = finalRes+ ω^-i*Pick[solFunctionsPower ,solFunctionsPower [[;;,2]],i] /. sol[[1]];
,{i,1,orderNo}]

Print["Simplifying..."];
results = Simplify[finalRes];

resList = Thread[solFunctions ==  Collect[ Limit[results,s-> Infinity],ω, FullSimplify]] /.s-> Infinity;

(*Assertion message definition*)
AbortAssert::trace="Non-zero solutions for higher Fourier components present.";
AbortAssert[test_]:=Check[TrueQ[test]||Message[AbortAssert::trace],Abort[]]

(* Check that at the end of the flow all non-zero Fourier components have vanished and print the results. 
   The obtained coefficients h_j multiply the generators G_j.
*)
assertFlag = True;
Do[If[i[[1]][[1]]!=0,If[i[[2]]=!=0,assertFlag = False; Break[]]],{i,resList}];
If[assertFlag == True, Print["Solutions:"]; Do[If[i[[1]][[1]]==0,Print[TraditionalForm[(Subscript[h, Position[resList,i][[1,1]]])==i[[2]]]]],{i,resList}],AbortAssert[assertFlag]];

(* Run this cell to obtaind the expansion for the micromotion, that is the operator S such that U = exp(-i S) as in Eq. (25).
The obtained coefficients d_j multiply the generators G_j.
*)
resArray = ArrayReshape[results,{fourierCompNo+1,basisNo}];
resArray = Expand[ϵ*resArray /. ω->1/ϵ]/.(ϵ^pow_/;pow>magnusTermNo)->0;
Do[fH[i]=resArray[[i+1]];,{i,0,fourierCompNo}]

Fix=Now;
A = Sum[(fH[n]*Exp[I*n*ω*t]-ComplexExpand[Conjugate[fH[n]],complexParameters]*Exp[-I*n*ω*t]),{n,1,fourierCompNo}];
Print["Now we have the generator ... "];

Print["doing 1 ... ", TimeObject[Fix]];
Print["Integrating ... ", Now-Fix];
Ω[1] = Integrate[A,{s,0,s1}]  /. s1->s;
magnusRes = Ω[1];

Do[
Print["doing ", n, " ... ", Now-Fix];
S[n,1] = Expand[Comm[Ω[n-1],A]] /. (ϵ^pow_/;pow>magnusTermNo)->0;
Do[
S[n,j] = Expand[Sum[Comm[Ω[m],S[n-m,j-1]],{m,1,n-j}]] /. (ϵ^pow_/;pow>magnusTermNo)->0;
,{j,2,n-1}];
Print["Integrating ... ", Now-Fix];
Ω[n] = Sum[BernoulliB[j]/j!*Integrate[S[n,j],{s,0,s1}],{j,1,n-1}] /.s1->s;
magnusRes = magnusRes +  Ω[n];
,{n,2,magnusTermNo}]

Print["proceeding to limits ... ", Now-Fix];
(*cutoff higher orders*)
magnusLim = Limit[magnusRes,s-> +Infinity];
Print["collecting ... ", Now-Fix];
resMagnusCol = Collect[magnusLim*(-I),ϵ,FullSimplify]/. ϵ -> 1/ω;
Do[Print[TraditionalForm[(Subscript[S, Position[resMagnusCol,i][[1,1]]]) == i]],{i,resMagnusCol}];
Print["All done ... ", Now-Fix];