Generating system of equations with unique solutions

If you are looking at intersecting the zero sets of two polynomials in two variables, then by Bezout's theorem you should anticipate that the number of solutions will be equal to the product of the degrees of the polynomials. So in general, you will not get unique solutions. You will need to be very careful to hide these extra solutions "at infinity" or in the complex plane, or as "double roots".

Generating systems. The same method that works for linear equations works also for polynomial equations. Starting with a solution in mind (in mathematics and computer science, we call this a planted solution), generate in some way left-hand sides of equations, and then compute the corresponding right-hand sides. If you use more than one equation per variable then you are likely not to get additional solutions.

Verifying uniqueness. Modern computer algebra systems can find the complete set of solutions of a system of polynomial equations. Under the hood, they might use Gröbner bases.

As an example, in sage you can write

sage: var('x, y')
(x, y)
sage: solve([x + y == 4, 2*x*y == 6], [x, y])
[[x == 1, y == 3], [x == 3, y == 1]]

Sage correctly found the two unique solutions.

Assuming that sage is using Gröbner bases, here is what happens. First, sage computes a Gröbner basis in elimination order.

sage: R.<x,y> = PolynomialRing(ZZ, 2, order = 'lex')
sage: I = ideal(x + y - 4, 2*x*y - 6)
sage: I.groebner_basis()
[x + y - 4, 2*y^2 - 8*y + 6]

It solves the quadratic to find all possible values of $y$, substitutes them, and finds the corresponding values of $x$. The same idea works with more variables, although the degrees can get quite large.