Is this type of nonlinear system known to factorize N ? What is the best way (computationally speaking) to solve it ?

edit

If You square eq1, then all your equations are multivariate quadratic polynomials, so you can compute a Gröbner basis to solve the system.

You can achieve this with the following code:

P.<N, x, y, x1, x2, x3, x4, y1, y2, y3> = PolynomialRing(QQ, order="degrevlex")

eq0 = (4*35+2)^2-(2*21-1)^2 -N
eq1 = 16*x^2+16*x+4-N-(2*y-1)^2

eq2 = (N+(2*y-1)^2-4)/4+1-2*(2*x+1) +1-(4*x1+2)^2
eq3 = ((4*x1+2)^2-4)/4+1-2*(2*x1+1) +1-(4*x2+2)^2
eq4 = ((4*x2+2)^2-4)/4+1 -2*(2*x2+1) +1-(4*x3+2)^2
eq5 = ((4*x3+2)^2-4)/4+1 -2*(2*x3+1) +1-(4*x4+2)^2
eq6 = ((4*x4+2)^2-4)/4+1-2*(2*x4+1) +1-36

eq7 = (N+(2*y-1)^2-16*(x-y)*(x+y+1)-4)/4+1-2*(2*y+1) +1-(4*y1+2)^2
eq8 = ((4*y1+2)^2-4)/4+1-2*(2*y1+1) +1-(4*y2+2)^2
eq9 = ((4*y2+2)^2-4)/4+1 -2*(2*y2+1) +1-(4*y3+2)^2
eq10 = ((4*y3+2)^2-4)/4+1-2*(2*y3+1) +1-36 

I = ideal([eq0, eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9, eq10])

gb = I.groebner_basis(algorithm="giac:gbasis", prot=True)

print(gb)
[N - 18483, x - 35, y - 21, x1 - 17, x2 + 9, x3 + 5, x4 + 3, y1 + 11, y2 - 5, y3 + 3]

For the complexity of solving the system, the example you provide consists of 9 quadratic polynomials in 10 variables. The DRL leading monomials of your input equations are:

for poly in I.gens():
    print(poly.lm())

N
x^2
y^2
x1^2
x2^2
x3^2
x4^2
x^2
y1^2
y2^2
y3^2

So your equations already contain a zero-dimensional DRL Gröbner basis, by ignoring one of the polynomials with leading monomial $x^2$. I am assuming here that this holds for the general case of $i + j + 1$ polynomials in $i + j$ variables too. You now have two options to solve the system:

1. Ignore one polynomial and perform term order conversion to LEX. With fast linear algebra this approach has polylogarithmic complexity (assuming that some precomputations are for free + theoretical assumptions I am not going to explain here) $\tilde{\mathcal{O}} \left( 2^{\omega \cdot (i + j)} \right)$.

2. You recompute the DRL Gröbner basis for the full system. With this and this paper the complexity of recomputing the DRL Gröbner basis is $\mathcal{O} \left( \binom{2 \cdot (i + j +1)}{i + j + 2}^\omega \right)$.

The constant $2 \leq \omega \leq 3$ is the exponent of matrix multiplication. State-of-the-art upper bound is $\omega < 2.371339$.