Output of CoefficientRules into Jacobian matrix

Question

I have a pair of large polynomials in u, v, a1...a9 and b1...b9, generated by

x1 := u
x2 := a1 u + a2 v + a3*u^2 + a4*u*v + a5*v^2 + a6 u^3 + a7 u^2 v + a8 u v^2 + a9 v^3 
y1 := v
y2 := b1*u + b2*v + b3*u^2 + b4*u*v + b5*v^2 + b6 u^3 + b7 u^2 v + b8 u v^2 + b9 v^3 
f1 := -(1 + k)*x1 - g*x1^3 + k*x2 - c x1 + c x2
f2 := k*x1 - k*x2 + c*x1 - c*x2

c := 0.2
k := 1
g := 0.7

poly1 = -y2 + a1*y1 + a2*f1 + 2*a3*x1*y1 + 2 a4*x1*f1 + a4*y1^2 + 2 a5*y1*f1 + 3 a6*x1^2*y1 + 2 a7*x1*y1^2 + a7 x1^2 f1 + a8 y1^3 +  2 a8 x1 y1 f1 + 3 a9 y1^2 f1
poly2 = -f2 + b1*y1 + b2*f1 + 2*b3*x1*y1 + 2 b4*x1*f1 + b4*y1^2 + 2 b5*y1*f1 + 3 b6*x1^2*y1 + 2 b7*x1*y1^2 + b7 x1^2 f1 + b8 y1^3 + 2 a8 x1 y1 f1 + 3 b9 y1^2 f1

I would like to be able to numerically approximate all the a* and b*. I need to set the coefficients of the u and v terms of up to order 3 to zero, then use something (Newton's method) to solve a system of 18 non-linear equations.

Using the following gives me the coefficients, but I don't know how to manipulate these results to give me what I want.

CoefficientRules[poly1, {u, v}]
CoefficientRules[poly2, {u,v}]

I now need to construct the [18]x[18] Jacobian and implement a non-linear solver.

Answer 1

eqs = Union @@ (Select[CoefficientRules[#, {u, v}], 
                       Tr@#[[1]] <= 3 &][[All, 2]] & /@ {poly1, poly2});
Reduce[Thread[eqs == 0]]

or

Solve[Thread[eqs == 0], Variables@eqs] // N