Consider a system $$ \sum_{0\le i,j\le N} x^i y^j A_{ij}(u,v) = 0 \qquad (*) $$ where $x,y,u,v\in \mathbb{R}$ and $A_{ij}$ are $3\times 3$ matrices, whose entries are smooth functions of wrt $(u,v)$.
Let $(u^*(x,y),v^*(x,y))$ be a solution of $(*)$, that may be hard to compute explicitly.
Is it possible to compute a Taylor expansion of $(u^*(x,y),v^*(x,y))$ at the point $(0,0)$ up to order $N$ wrt $x,y$ ?
Consider the case $N=1$, so that the approximation we would like to compute is a plane. Rewrite $(*)$ as $F(x,y,u,v)=0$. Applying the IFT, we can find a neighborhood $N$ of $0$ and a function $g$ such that $F(x,y,g(x,y))=0$ in $N$ and $$ \frac{\partial g}{\partial u}(0,0) = -Jac(F)_{xy}(x,y,0,0)^{-1}\ \frac{\partial F}{\partial u}(x,y,0,0) \qquad (1) \\ \frac{\partial g}{\partial v}(0,0) = -Jac(F)_{xy}(x,y,0,0)^{-1}\ \frac{\partial F}{\partial v}(x,y,0,0)\qquad (2) $$ and then the approximation would be given by $$ (u^*,v^*)^\intercal = g(0,0) + u \frac{\partial g}{\partial u}(0,0) + v \frac{\partial g}{\partial v}(0,0) + o(u,v) $$