Program por computing the relevant expressions related with the Wilcox expansion such as they appear in the paper “Sequences of transformations and perturbative exponential factorizations”, by A. Arnal, F. Casas, C. Chiralt and J.A. Oteo
Date: 30/10/2020
To run the program, just Evaluate Notebook in the menu Evaluation.  By default, the program computes only the first 7 orders (in the Results), but higher orders are obtained in a trivial way.

Initialization

This block defines the properties of the commutator (Cmt)

$RecursionLimit=Infinity;

$IterationLimit=Infinity;

Cmt, ad

Cmt/:Format[Cmt[a_,b_]]:=SequenceForm["[",a,",",b,"]"];

ad[a_,0,b_]:=b;
ad[a_,j_Integer,b_]:=Cmt[a,ad[a,j-1,b]];

Cmt[x_,x_]:=0;
Cmt[a___, 0,     b___] := 0;
Cmt[a___, x_+y_, b___] :=    Cmt[a,x,b]+Cmt[a,y,b];
Cmt[a___, t^n_ x_,  b___] := t^n Cmt[a,x,b];
Cmt[a___, n_ x_Cmt, b___] := n   Cmt[a,x,b];
Cmt[a___, x_  y_,  b___] :=  Cmt[a, Expand[x y], b];

Next we define a non-commutative product (Mm) by assigning the required properties

Mm[a___, x_+y_, b___] :=    Mm[a,x,b]+Mm[a,y,b];
Mm[a___, h^n_ x_,  b___] := h^n Mm[a,x,b];
Mm[a___, x_ y_Mm, b___] := x   Mm[a,y,b];
Mm[a___, x_  y_,  b___] :=  Mm[a, Expand[x y], b];
Mm[a___, n_ ide, b___] := n * Mm[a,ide,b];
SetAttributes[Mm,Flat];

Mm[{a___}]:=Mm[a]

Expressions for (Subscript[W, n])(t)

Recursive procedure, Eqs. (11)-(16): (Subscript[W, n]) in terms of (Subscript[W, 1, ])...(Subscript[W, n-1 ])

b[0,1]:=A[t];b[0,n_]:=0;
b[n_,l_]:=b[n,l]=Sum[((-1)^k/k!)ad[W[n],k,b[n-1,l-n k]-g[n,l-n k]],{k,0,IntegerPart[(l-1)/n]-1,1}]//Simplify;
g[n_,m_]:=0;g[1,k_]:=g[1,k]=(1/k!)ad[W[1],k-1,b[0,1]];
g[n_,r_]:=(1/(r/n)!)ad[W[n],(r/n)-1,b[n-1,n]]/;IntegerQ[r/n];
Wdot[n_]:=Wdot[n]=Expand[b[n-1,n]];

Results up to n=7:

Wdot[1]

Wdot[2]

Wdot[3]

Wdot[4]

Wdot[5]

Wdot[6]

Wdot[7]

Explicit expresion of (Subscript[P, n]) in terms of (Subscript[W, 1])...(Subscript[W, n ]) Equation (19)

P[n_]:=W[n]+Sum[Mm[Map[W,Reverse[IntegerPartitions[n][[s]] ]]]/(Product[Counts[IntegerPartitions[n][[s]]][[j]]!,{j,1,Length[Counts[IntegerPartitions[n][[s]]]]}]),{s,2,Length[IntegerPartitions[n]]}]

Results up to n = 7:

P[1]

P[2]

P[3]

P[4]

P[5]

P[6]

P[7]

By working out Equation (19) one can invert the relations and express (Subscript[W, n ]) in terms of (Subscript[P, 1])...(Subscript[P, n])

Clear[P,W]

W[1]:=P[1]

W[n_]:=Mm[P[n]]-Sum[Mm[Map[W,Reverse[IntegerPartitions[n][[s]] ]]]/(Product[Counts[IntegerPartitions[n][[s]]][[j]]!,{j,1,Length[Counts[IntegerPartitions[n][[s]]]]}]),{s,2,Length[IntegerPartitions[n]]}]//Simplify//Expand

Results up to n = 7 :

W[2]

W[3]

W[4]

W[5]

W[6]

W[7]

Now after doing these (Subscript[P, i]) products we get Equation (28), that is,(Subscript[W, n ])in terms of independent nested commutators of A(t)

Code for doing the product

ww[n_]:=ww[n]=W[n]/.Mm->h

h[a___,b_,c_]:=h[a,dott[b,c]]

Now we compute and rename the first ten (Subscript[W, n ]). More terms must be computed and renamed if one wants to compute more (Subscript[W, n ]) in terms of nested independent commutators of A(t).

Do[w[n]=ww[n],{n,1,10}];

dott[string1_,string2_]:=dott[string1,string2]=
If[Length[string1[[1]]]==0,If[Length[string2[[1]]]==0,dot[string1,string2],pdot[string1,string2]],If[Length[string2[[1]]]==0,Union[Flatten[Table[dot[string1[[i]],string2],{i,1,Length[string1]}],1 ]],Union[Flatten[Table[pdot[string1[[i]],string2],{i,1,Length[string1]}],1 ]]]]

pdot[string1_,string2_]:=pdot[string1,string2]=Union[Flatten[Table[dot[string1,string2[[i]]],{i,1,Length[string2]}],1 ]]

dot[string1_,string2_]:=dot[string1,string2]=Module[{p,tp,ord,p2,s1,s2,shuf},s1=InversePermutation[string1];s2=InversePermutation[string2];p2=s2+Length[s1];p=Permutations@Join[1&/@s1,0&/@p2]\[Transpose];
tp=BitXor[p,1];
ord=Accumulate[p] p+(Accumulate[tp]+Length[s1]) tp;
shuf=Outer[Part,{Join[s1,p2]},ord,1][[1]]\[Transpose];
Table[InversePermutation[shuf[[j]]],{j,1,Length[shuf]}]]

P[n_]:=Table[i,{i,1,n}]

h[list_]:=If[Length[list[[1]]]==0,f[list],Sum[Apply[f,{list[[i]]}],{i,1,Length[list]}]]

Equation (28), with A((Subscript[t, 1])) as the last operator

f[{list___}]:=0;

f[{list___,1}]:=A[{list,1}];

W1[n_]:=W1[n]=w[n]//Simplify//Expand

Results up to n = 5 :

W1[2]

W1[3]

W1[4]

W1[5]

Equation (28), with A((Subscript[t, n])) as the last operator

Clear[f]

f[{list___}]:=f[Length[{list}]][{list}]

f[x_][{list___}]:=0;

f[x_][{list___,x_}]:=A[{list,x}];

Wn[n_]:=Wn[n]=w[n]//Simplify//Expand

Results up to n = 5 :

Wn[2]

Wn[3]

Wn[4]

Wn[5]

Notice that the number of non-vanishing terms is different in each representation