# Best way to determine polynomial coefficients in series expansion

I would like to solve the equation $$h'(\boldsymbol{x}_1)\left[B_1\boldsymbol{x}_1+g_1(\boldsymbol{x_1},h(\boldsymbol{x}_1))\right]=B_2h(\boldsymbol{x}_1)+g_2(\boldsymbol{x}_1,h(\boldsymbol{x}_1))$$

where $B_1$ is a matrix of $\mathbb{R}^2$, $B_2$ is a "matrix" of $\mathbb{R}$, $g_1:\mathbb{R}^3\longrightarrow \mathbb{R}^2$, $g_2:\mathbb{R}^3\longrightarrow\mathbb{R}$.

$h$ is the unknow function, search in terms of its Taylor expansion in a neighbourhood of 0: $$h(u,v)\approx\sum_{k=0}^d \sum_{l=0}^k a_{kl}u^k\,v^l$$

I injected the second equation in the left-hand side of the first equation, denoted the result by left, and injected again the second equation in the right-hand side to obtain right.

Given left and right, I want to calculate the $a_{ij}$, when it is possible (in $\mathbb{R}$).

I tried several solutions:

1. use CoefficientList to extract the polynomial coefficients of right-left and solve ==0 for each term.

2. take every derivative in $u$ and $v$ and take the values in $0$, which is probably the same as what CoefficientList does, but manually.

3. evaluate left and right for many different values of $u$ and $v$ and solve the system.

Which solution would you use, and why? In particular, how would you treat the terms of degre >$d$ when calculating the $a_{ij}$?

-