Inverting a system of series

The problem is as follows, I have some series of expansion $$t_1=\sum_{n_1,n_2}a^{(1)}_{n_1,n_2}x_1^{n_1}x_2^{n_2},\quad t_2=\sum_{n_1,n_2}a^{(2)}_{n_1,n_2}x_1^{n_1}x_2^{n_2}$$ and I want to express the variables $$x_1$$ and $$x_2$$ in terms of the result of the sums $$t_1,t_2$$. i.e I want to find $$b^{(1)}_{n_1,n_2},b^{(2)}_{n_1,n_2}$$ such that, $$x_1=\sum_{n_1,n_2}b^{(1)}_{n_1,n_2}t_1^{n_1}t_2^{n_2},\quad x_2=\sum_{n_1,n_2}b^{(2)}_{n_1,n_2}t_1^{n_1}t_2^{n_2}$$.

If I had just one variable $$t=\sum a_n x^n$$ InverseSeries will do the trick, but having several variables seems to make the problem trickier.

There are a couple of remarks to have into account, first, this is of course not possible for arbitrary series, but the series I have start as $$t_i=x_i+\text{higher order}$$ and in the particular problem at hand I am sure the inversion is possible. Second, the two variables are of the same order, i.e I want to take all terms in the series with $$n_1+n_2\leq N$$.

Now, what I have tried (I think is not a very good solution, so feel free to ignore it): As mathematica has some tricky behaviour for several parameters series (it seems like the only way to get is to use nested series), I introduce some parameter $$\epsilon$$ allowing me to control the expansion and set $$t_i\to \epsilon t_i,x_i\to\epsilon x_i$$. Then, already truncating, the series are $$t_i\epsilon=\sum_{0\leq n_1,n_2\leq N}a^{(i)}_{n_1,n_2}(\epsilon x_1)^{n_1}(\epsilon x_2)^{n_2}+O[\epsilon]^{N+1}$$ Then, the only thing that occur to me, was to make the expansion of $$x_i$$ as $$x_i\epsilon=\sum_{0\leq n_1,n_2\leq N}b^{(i)}_{n_1,n_2}(\epsilon t_1)^{n_1}(\epsilon t_2)^{n_2}+O[\epsilon]^{N+1}$$ and plug it in in the original series, Then use LogicalExpand to equate the resulting series to $$\epsilon t_i$$. This gives me some equation for the $$b$$'s for which I can then use Solve.

There are two things that I don't like about this solution, first of all, it is too much work: I have been working with 1d systems (i.e just one $$t$$ and one $$x$$) and InverseSeries works like a charm, now in 2d systems it seems like an easy calculation became horribly complicated and I will need to lower considerably the order of the series ($$N$$) to get results in some decent time. Moreover, I would like to go to higher dimensional systems ($$t_1,\dots,t_n$$ and $$x_1,\dots,x_n$$) and it seems that things will get incredibly messy. Second, The logical expansion is done in this auxiliary parameter $$\epsilon$$ and not in monomials $$t_1^{n_1}t_2^{n_2}$$; this not only looks like a terrible practice, but it also seems very likely that Mathematica at some point will just give me some solution of $$b$$ as functions of t.

Here the program, I am sorry to change the notation, but I have been puzzling for a while now and I was running out of letters, hope the explanation above was clear. (Note that in fact it doesn't give a good result because it solves for $$d1,d2$$ that are not to be solved). There you will also see the explicit values for the series at $$N=3$$.

In[83]:= u1 = Exp[t1ser]

Out[83]= SeriesData[s, 0, {
x1, x1 (92 x1 - x2),
Rational[1, 2]
x1 ((92 x1 - x2)^2 + 2 (
3812 x1^2 + 46 x1 x2 + Rational[-3, 2] x2^2))}, 1, 4, 1]

In[84]:= u2 = Exp[t2ser]

Out[84]= SeriesData[s, 0, {
x2, x2 (46 x1 + 2 x2),
Rational[1, 2]
x2 ((46 x1 + 2 x2)^2 + 2 (
Rational[9599, 2] x1^2 - 92 x1 x2 + 3 x2^2))}, 1, 4, 1]

In[95]:= try = {x1 ->
Sum[a[n1, n2] (s d1)^n1 (s d2)^n2, {n1, 0, M}, {n2, 0, M - n1}] ,
x2 -> Sum[
b[n1, n2] (s d1)^n1 (s d2)^n2, {n1, 0, M}, {n2, 0, M - n1}]}

Out[95]= {x1 ->
a[0, 0] + d2 s a[0, 1] + d2^2 s^2 a[0, 2] + d1 s a[1, 0] +
d1 d2 s^2 a[1, 1] + d1^2 s^2 a[2, 0],
x2 -> b[0, 0] + d2 s b[0, 1] + d2^2 s^2 b[0, 2] + d1 s b[1, 0] +
d1 d2 s^2 b[1, 1] + d1^2 s^2 b[2, 0]}

In[81]:= var =
Flatten[Table[{a[n1, n2], b[n1, n2]}, {n1, 0, M}, {n2, 0, M - n1}]]

Out[81]= {a[0, 0], b[0, 0], a[0, 1], b[0, 1], a[0, 2], b[0, 2],
a[1, 0], b[1, 0], a[1, 1], b[1, 1], a[2, 0], b[2, 0]}

In[120]:= coef1 = LogicalExpand[d1 s == u1 /. try];
coef2 = LogicalExpand[d2 s == u2 /. try];

In[122]:= coef = Join[coef1, coef2]

Out[122]= -d1 + a[0, 0] == 0 &&
92 a[0, 0]^2 + d2 a[0, 1] + d1 a[1, 0] - a[0, 0] b[0, 0] == 0 &&
8044 a[0, 0]^3 + 184 d2 a[0, 0] a[0, 1] + d2^2 a[0, 2] +
184 d1 a[0, 0] a[1, 0] + d1 d2 a[1, 1] + d1^2 a[2, 0] -
46 a[0, 0]^2 b[0, 0] - d2 a[0, 1] b[0, 0] - d1 a[1, 0] b[0, 0] -
a[0, 0] b[0, 0]^2 - d2 a[0, 0] b[0, 1] - d1 a[0, 0] b[1, 0] ==
0 && -d2 + b[0, 0] == 0 &&
46 a[0, 0] b[0, 0] + 2 b[0, 0]^2 + d2 b[0, 1] + d1 b[1, 0] == 0 &&
11715/2 a[0, 0]^2 b[0, 0] + 46 d2 a[0, 1] b[0, 0] +
46 d1 a[1, 0] b[0, 0] + 5 b[0, 0]^3 + 46 d2 a[0, 0] b[0, 1] +
4 d2 b[0, 0] b[0, 1] + d2^2 b[0, 2] + 46 d1 a[0, 0] b[1, 0] +
4 d1 b[0, 0] b[1, 0] + d1 d2 b[1, 1] + d1^2 b[2, 0] == 0

In[123]:= Solve[coef]

Out[123]= {{d1 -> 0, d2 -> 0, a[0, 0] -> 0,
b[0, 0] -> 0}, {d1 -> a[0, 0], d2 -> b[0, 0],
a[1, 0] -> (-92 a[0, 0]^2 + a[0, 0] b[0, 0] - a[0, 1] b[0, 0])/
a[0, 0],
a[2, 0] -> (1/(
a[0, 0]^2))(8884 a[0, 0]^3 - 276 a[0, 0]^2 b[0, 0] -
a[0, 0] a[1, 1] b[0, 0] - a[0, 2] b[0, 0]^2),
b[1, 0] -> -((b[0, 0] (46 a[0, 0] + 2 b[0, 0] + b[0, 1]))/a[0, 0]),
b[2, 0] -> (1/(2 a[0, 0]^2))
b[0, 0] (981 a[0, 0]^2 + 460 a[0, 0] b[0, 0] + 6 b[0, 0]^2 -
2 b[0, 0] b[0, 2] - 2 a[0, 0] b[1, 1])}, {d1 -> 0,
d2 -> b[0, 0], a[0, 0] -> 0, a[0, 1] -> 0, a[0, 2] -> 0,
b[0, 1] -> -2 b[0, 0], b[0, 2] -> 3 b[0, 0]}}


Perhaps you can use the new in M12 function AsymptoticSolve. For example:

eq = And[
t1 == Sum[b1[i,j] x1^i x2^j, {i, 0, 2}, {j, 0, 2}] - b1[0,0],
t2 == Sum[b2[i,j] x1^i x2^j, {i, 0, 2}, {j, 0, 2}] - b2[0,0]
]


t1 == x2 b1[0, 1] + x2^2 b1[0, 2] + x1 b1[1, 0] + x1 x2 b1[1, 1] + x1 x2^2 b1[1, 2] + x1^2 b1[2, 0] + x1^2 x2 b1[2, 1] + x1^2 x2^2 b1[2, 2] && t2 == x2 b2[0, 1] + x2^2 b2[0, 2] + x1 b2[1, 0] + x1 x2 b2[1, 1] + x1 x2^2 b2[1, 2] + x1^2 b2[2, 0] + x1^2 x2 b2[2, 1] + x1^2 x2^2 b2[2, 2]

Then:

AsymptoticSolve[eq, {{x1, x2}, {0,0}}, {{t1, t2}, {0,0}, 1}] //TeXForm


$$\left\{\left\{\text{x1}\to \frac{\text{t1} \text{b2}(0,1)}{\text{b1}(1,0) \text{b2}(0,1)-\text{b1}(0,1) \text{b2}(1,0)}+\frac{\text{t2} \text{b1}(0,1)}{\text{b1}(0,1) \text{b2}(1,0)-\text{b1}(1,0) \text{b2}(0,1)},\text{x2}\to \frac{\text{t1} \text{b2}(1,0)}{\text{b1}(0,1) \text{b2}(1,0)-\text{b1}(1,0) \text{b2}(0,1)}+\frac{\text{t2} \text{b1}(1,0)}{\text{b1}(1,0) \text{b2}(0,1)-\text{b1}(0,1) \text{b2}(1,0)}\right\}\right\}$$

I think the system of equations you use in your example is:

eqs = And[
t1 == x1+92 x1^2+8044 x1^3-x1 x2-46 x1^2 x2-x1 x2^2,
t2 == x2+46 x1 x2+(11715 x1^2 x2)/2+2 x2^2+5 x2^3
];


and the solution returned by AsymptoticSolve is:

AsymptoticSolve[eqs, {{x1, x2}, {0, 0}}, {{t1, t2}, {0, 0}, 3}] //TeXForm


$$\left\{\left\{\text{x1}\to \text{t1}-92 \text{t1}^2+8884 \text{t1}^3+\text{t1} \text{t2}-276 \text{t1}^2 \text{t2},\text{x2}\to \text{t2}-46 \text{t1} \text{t2}+\frac{981 \text{t1}^2 \text{t2}}{2}-2 \text{t2}^2+230 \text{t1} \text{t2}^2+3 \text{t2}^3\right\}\right\}$$