I have some large system of particular non-linear polynomial equations:

- each equation
*mentions at most three variables* - no variable appears with a degree larger than 1.

I'm not an expert in this area but I think these are called "multilinear" polynomials, and in particular "trilinear" in my case, but I have the additional constraint that there are only three variable mentioned in each equation.

For example:

\begin{equation} \begin{aligned} 3x_1x_2x_3 + 2x_2x_3 - 3x_1x_3 + 1 &= 0\\ \ldots &= \ldots \\ 2x_1x_{42}x_{57} - 3x_1x_{57} - 2x_{42} + 2 & = 0 \end{aligned} \end{equation}

Can we "efficiently" solve these systems numerically? I know non-linear systems are hard but I've read multilinear polynomials have a lot of special properties so I suppose solving methods can benefit from ad-hoc techniques. And then, I have this additional constraint on the number of variables mentioned in each equation, which may help, if at all, at least to represent the system sparsely.

Can someone point me to relevant literature, if any, on these kind of systems with particular focus on numerical methods for solving them?

**Edit**: from the comments I realized the question could be more specific. Of course the obvious answer is Newton's method. However what I would like to understand is whether the particular shape of my equations can:

- help the method converge faster, or
- help the method be faster in practice/be more efficiently implementable (e.g. because of a better way to store the coefficients), or
- enable some different method entirely

Another detail: I expect the system to underconstrained most of the times.