We have a system of $m$ bilinear equations in $\mathrm x, \mathrm y \in \mathbb R^n$

$$\begin{aligned} \mathrm x^\top \mathrm A_1 \,\mathrm y &= b_1\\ \mathrm x^\top \mathrm A_2 \,\mathrm y &= b_2\\ &\vdots\\ \mathrm x^\top \mathrm A_m \,\mathrm y &= b_m\end{aligned}$$

If matrices $\mathrm A_1, \mathrm A_2, \dots, \mathrm A_m$ are very sparse, perhaps it would not be utterly hopeless to use symbolic methods (e.g., Gröbner bases). However, we can use numerical methods. Note that

$$b_i = \mathrm x^\top \mathrm A_i \,\mathrm y = \mbox{tr} \left( \mathrm x^\top \mathrm A_i \,\mathrm y \right) = \mbox{tr} \left( \mathrm y^\top \mathrm A_i^\top \,\mathrm x \right) = \mbox{tr} \left( \mathrm A_i^\top \mathrm x \mathrm y^\top \right) =: \langle \mathrm A_i, \mathrm x \mathrm y^\top\rangle$$

Let $\mathrm Z := \mathrm x \mathrm y^\top$. Hence, we have $m$ linear equality constraints in $\rm Z$ and a constraint on its rank

$$\begin{aligned} \langle \mathrm A_1, \mathrm Z \rangle &= b_1\\ \langle \mathrm A_2, \mathrm Z\rangle &= b_2\\ &\vdots\\ \langle \mathrm A_m, \mathrm Z\rangle &= b_m\\ \mbox{rank} (\mathrm Z) &= 1\end{aligned}$$

Since the nuclear norm is a convex proxy for the rank, we solve the following convex program in $\rm Z$

$$\begin{array}{ll} \text{minimize} & \| \mathrm Z \|_*\\ \text{subject to} & \langle \mathrm A_1, \mathrm Z \rangle = b_1\\ & \langle \mathrm A_2, \mathrm Z \rangle = b_2\\ & \qquad\vdots\\ & \langle \mathrm A_m, \mathrm Z \rangle = b_m\end{array}$$

If the optimal solution is rank-$1$, then we have solved the original system of $m$ bilinear equations.