The problem of inverting the implicit function $$y=f(x)$$ in the form of power-series $$x = a + \sum_{k=1}^{\infty} b_k (y-f(a))^k$$
around a point of interest $x=a$, has a long history. Lagrange obtained a theoretical inversion formula ((https://mathworld.wolfram.com/LagrangeInversionTheorem.html)), yet efficient implementations are relatively recent ((Brent, Richard P., and Hsiang T. Kung. “Fast algorithms for manipulating formal power series.” Journal of the ACM (JACM) 25, no. 4 (1978): 581-595.)) ((Johansson, F., 2015. A fast algorithm for reversion of power series. Mathematics of Computation, 84(291), pp.475-484.)).
I this note I am sharing a bit simpler algorithm, performing Newton-like updates:
$$x_{k} = x_k-f(x_{k-1})[(y-f(a))^{k}]\cdot \frac{(y-f(a))^{k}}{f'(a)},\quad x_1 = a+\frac{y-f(a)}{f'(a)}.$$
We can see that it gradually produces more and more accurate terms. More precisely, assuming w.l.o.g. \(0=a=f(a)\), suppose \(f(x) =y + O(y^{k})\), then \(f(x + h ) = f(x)+f'(0) h + O(h^2)\) by Taylor’s expansion, and so the update \(h=-\frac{f(x) [y^k] }{f'(0)}\) (where \([z^k]\) is the operation of extracting the coefficient with \(z^k\)) gives \(f(x+h)=y+O(y^{k+1})\) as desired.
I implemented this a proposal for Sympy.