I'm looking for an algorithm that can quickly generate an approximate real number solutions of an underdetermined quadratic system. My task is to explore its algebraic variety. I'm interested in a special cases of quadratic systems that arise from mechanics:

set of variables $X = \{ x\_i, y\_i, z\_i | 1 \leqslant i \leqslant N\}$ (3\*N variables for $N \in \mathbb{N}$),

a system consists of two kinds of equations:

1) $(x\_i-x\_j)^2+(y\_i-y\_j)^2+(z\_i-z\_j)^2 + p\_{ij} = 0$ for $i, j \in [1,N]$ and $p\_{ij} \in \mathbb{R}$.

2) a more general case of 1: $(x\_i-x\_k)(x\_j-x\_k)+(y\_i-y\_k)(y\_j-y\_k)+(z\_i-z\_k)(z\_j-z\_k) + a\_{ijk} = 0 $ for $i, j, k \in [1,N]$ and $a\_{ijk} \in \mathbb{R}$.

It's guaranteed that variety of solutions is stable under small changes of coefficients.
Here is an example with six-dimensional variety of solutions:

```
(x1-x2)^2 + (y1-y2)^2 + (z1-z2)^2 - 1 = 0
(x3-x2)^2 + (y3-y2)^2 + (z3-z2)^2 - 1 = 0
(x3-x4)^2 + (y3-y4)^2 + (z3-z4)^2 - 1 = 0
(x1-x4)^2 + (y1-y4)^2 + (z1-z4)^2 - 1 = 0
(x2-x1)*(x4-x1) + (y2-y1)*(y4-y1) + (z2-z1)*(z4-z1) - 1/3 = 0
(x1-x2)*(x3-x2) + (y1-y2)*(y3-y2) + (z1-z2)*(z3-z2) - 1/3 = 0
(x2-x3)*(x4-x3) + (y2-y3)*(y4-y3) + (z2-z3)*(z4-z3) - 1/3 = 0
(x1-x4)*(x3-x4) + (y1-y4)*(y3-y4) + (z1-z4)*(z3-z4) - 1/3 = 0

```

First, the algorithm is allowed to process this underdetermined system for some time. After that, it should get as input a set of constraints, like $x\_1 = 0$ or $z\_N = 1$, that turn the system into determined and quickly generate real solutions of it. For instance, one can find 8 real solutions of the system above after constraining all of $\{x\_4,y\_4,z\_4,y\_3,z\_3,z\_2\}$ to zero. 
**I need the algorithm to be scalable** so it can be applied for N=25-30, however **high precision of solutions isn't required**, 3-4 orders will do.

I've come up with a few approaches:

1. I tried to calculate Gröbner bases for such systems in Magma:

```
    P< x1,y1,z1,... > := PolynomialRing(RationalField(),3*N);
    I := ideal<P | my equations >;
    Groebner(I);

```

It works for zero-dimensional ideals, but coefficients were big even for N=7. In cases with dimensions >0, calculation took too long (with enough RAM) and I have been able to compute Gröbner bases only for small systems.
So I tried compute approximate Gröbner bases by running OpenF4 algorithm with datatype of coefficients changed on double, but the error grew significantly on each step of the computation.

2. Also I consider a calculation of Gröbner basis in grevlex order with subsequent transformation to lexicographic using FGLM-like algorithm.
3. I found this article <https://arxiv.org/abs/1703.07810> and this looks useful for my purpose. I'm not familiar with such methods, but probably this is exactly what I need.

Here are my questions:

Is there a way to calculate an approximate Gröbner bases in reasonable time using my idea or any other approach, maybe for other order with subsequent transformation to lexicographic? 

Can this problem be solved without Gröbner bases using methods like №3 or any other?

Any comments on my ideas and new suggestions would be highly appreciated.