Help implementing Magnus Expansion

Question

The Magnus expansion is a tool to approximate solutions to first-order linear differential equations (the Wikipedia page is quite instructive and concise) - it's particularly useful because all orders in the expansion can be generated via recurrence relations.

However, a problem arises because as part of this recurrence you need to perform a series of higher-order integrals with indexed integration variables - I'm having a hard time implementing a programmatic way to generate the Magnus expansion for arbitrary order because I can't think of a way to generate a recurrence relation with a nested integral whose integration bounds are variable (see for instance $ \Omega_4 $ on the Wikipedia page).

I'll attach my code below, but I'm not sure how instructive it will be, but any help at all on solving this problem would be greatly appreciated. In my research I did find a Mathematica package a few years out of date which seems to use explicitly subscripted variables to keep track of the nested indices - however I would prefer a "home-built" solution which I can really understand rather than using that older package.

My current code (note that I have two implementations of the actual magnus Omega function, neither of which work as I want):

MagnusInt[ks_, H_] := Integrate[Fold[FComm, H[τ], MagnusΩ[#, H] & /@ ks], {τ, 0, t}]

MagnusS[n_, j_, H_] := Which[
  j == 1, Comm[MagnusΩ[n - 1, H], H[τ]],
  j == n - 1, AdComm[MagnusΩ[1, H], H[τ]],
  True, Sum[
   Comm[MagnusΩ[m, H], MagnusS[n - m, j - 1, H]],
   {m, 1, n - j}]
  ]

MagnusΩ1[n_, H_] := If[n == 1,
  Integrate[H[τ], {τ, 0, t}],
  Sum[BernoulliB[j]/
    j! Integrate[MagnusS[n, j, H[τ]], {τ, 0, t}], {j, 1, 
    n - 1}]
  ]

MagnusKs[n_, j_] := 
 Flatten[Permutations[#] & /@ IntegerPartitions[n - 1, {j}], 1]

MagnusΩ2[n_, H_] := If[n == 1,
  Integrate[H[τ], τ],
  Sum[BernoulliB[j]/j!*
    Total[
     MagnusInt[#, H] & /@ MagnusKs[n, j]
     ]
   , {j, 1, n - 1}]
  ]

Comm[A_, B_] := A.B - B.A;
FComm[A_, B_] := B.A - A.B;
AdComm[A_, B_, k_] := If[k == 0, B, Comm[A, AdComm[A, B, k - 1]]]

A useful test case I have been working with is for the matrix X = {{2, t}, {0, -1}}, for which the first few orders of the magnus expansion are calculated and shown on Page 18 in this PDF.

Any direct code suggestions or guidance on directions to explore would be greatly appreciated. Thanks!