How to isolate the real roots of a hard multivariate polynomial system with a finite number of solutions completely?

Question

In this marketing literature, the Wolfram Research alleges that

Mathematica automates root finding behind the higher-level functions Reduce and NSolve, removing the need for low-level operations like Maple’s Isolate.

However, in certain cases, it is sufficient to find (disjoint) isolating intervals for ALL real roots of a polynomial system. In other words, there is no need for exact solutions given by the advanced command Solve or Reduce. For a univariate polynomial with arbitrary real coefficients, the RootIntervals is of course useful, and in accordance with the above official claim, there appears to be a built-in function that performs analogous functionality on multivariate polynomial systems with finitely many real solutions, but I still cannot find out such an internal function in Mathematica as yet.
One may believe that NSolveValues should be enough here. Unfortunately, this is not always guaranteed, since the returned roots can be incomplete (even of a real polynomial system) without any warning messages. (Accordingly, the answer to this old question is simply "no", while space does not permit the citation of such examples.) Note that although a complete list of symbolic solutions may be found, the elapsed time would be much long (in contrast to solving numerically), which is more or less unacceptable. So, how do I obtain all rational isolating intervals (or boxes) for the multivariate case over the reals in Mma efficiently?

Edit. Here are some instances.
The first polynomial system is from $eq. (14)$ with $l=3$ (alternatively, with $l=2$) of Generalized Virasoro constructions for SU(3) (page 589):

eqs1 = {-Subscript[x,1]+4*l*Subscript[x,1]^2+8*Subscript[x,1]*(Subscript[x,2]+Subscript[x,3])-8*Subscript[x,2]*Subscript[x,3]+2*Subscript[x,1]*(Subscript[x,4]+Subscript[x,5]+Subscript[x,6]+Subscript[x,7])-2*Subscript[x,4]*Subscript[x,7]-2*Subscript[x,5]*Subscript[x,6],
-Subscript[x,2]+4*l*Subscript[x,2]^2+8*Subscript[x,2]*(Subscript[x,1]+Subscript[x,3])-8*Subscript[x,1]*Subscript[x,3]+2*Subscript[x,2]*(Subscript[x,4]+Subscript[x,5]+Subscript[x,6]+Subscript[x,7])-2*Subscript[x,4]*Subscript[x,6]-2*Subscript[x,5]*Subscript[x,7],
-Subscript[x,3]+4*l*Subscript[x,3]^2+8*Subscript[x,3]*(Subscript[x,1]+Subscript[x,2])-8*Subscript[x,1]*Subscript[x,2]+2*Subscript[x,3]*(Subscript[x,4]+Subscript[x,5]+Subscript[x,6]+Subscript[x,7])-2*Subscript[x,4]*Subscript[x,5]-2*Subscript[x,6]*Subscript[x,7],
-Subscript[x,4]+4*l*Subscript[x,4]^2+2*Subscript[x,1]*(Subscript[x,4]-Subscript[x,7])+2*Subscript[x,2]*(Subscript[x,4]-Subscript[x,6])+2*(Subscript[x,3]+3*Subscript[x,8])*(Subscript[x,4]-Subscript[x,5])+2*Subscript[x,4]*(4*Subscript[x,5]+Subscript[x,6]+Subscript[x,7]),
-Subscript[x,5]+4*l*Subscript[x,5]^2+2*Subscript[x,1]*(Subscript[x,5]-Subscript[x,6])+2*Subscript[x,2]*(Subscript[x,5]-Subscript[x,7])-2*(Subscript[x,3]+3*Subscript[x,8])*(Subscript[x,4]-Subscript[x,5])+2*Subscript[x,5]*(4*Subscript[x,4]+Subscript[x,6]+Subscript[x,7]),
-Subscript[x,6]+4*l*Subscript[x,6]^2-2*Subscript[x,1]*(Subscript[x,5]-Subscript[x,6])-2*Subscript[x,2]*(Subscript[x,4]-Subscript[x,6])+2*(Subscript[x,3]+3*Subscript[x,8])*(Subscript[x,6]-Subscript[x,7])+2*Subscript[x,6]*(Subscript[x,4]+Subscript[x,5]+4*Subscript[x,7]),
-Subscript[x,7]+4*l*Subscript[x,7]^2-2*Subscript[x,1]*(Subscript[x,4]-Subscript[x,7])-2*Subscript[x,2]*(Subscript[x,5]-Subscript[x,7])-2*(Subscript[x,3]+3*Subscript[x,8])*(Subscript[x,6]-Subscript[x,7])+2*Subscript[x,7]*(Subscript[x,4]+Subscript[x,5]+4*Subscript[x,6]),
-Subscript[x,8]+4*l*Subscript[x,8]^2+6*(Subscript[x,4]+Subscript[x,5]+Subscript[x,6]+Subscript[x,7])*Subscript[x,8]-6*Subscript[x,4]*Subscript[x,5]-6*Subscript[x,6]*Subscript[x,7],
(Subscript[x,3]+Subscript[x,8])*(Subscript[x,4]+Subscript[x,5]-Subscript[x,6]-Subscript[x,7])-2*Subscript[x,4]*Subscript[x,5]+2*Subscript[x,6]*Subscript[x,7]}/. l->3;

The second polynomial system is from $Example : CAPRASSE, DEMARET$ of Symbolic solution polynomial equation systems with symmetry (page 8):

eqs2=Array[3*Subscript[x,#]^3-Subscript[x,#]*(3*Subscript[x,#]*Indexed[a,1]+3*Indexed[a,1]+3*Indexed[a,2]-6)+Indexed[a,2]-Indexed[a,3]&,5]/.a:>Array[PowerSymmetricPolynomial[#,Array[Subscript[x,#]&,5]]&,3];(*$Assumptions:=Array[Subscript[x,#]&,8]\[Element]Reals*)

As for the isolation, these two instances ought to, ideally, take about sixteen seconds ….

Answer 1

There is one problem only, I see, in the second example: Subscripts.

      qs2 = qs1 /. {Subscript[a_, b_] :> 
              ToExpression[ToString[a] <> ToString[b]]};
  vars=Union@Flatten@Cases[qs2,_Symbol,Infinity]

     {x1, x2, x3, x4, x5, x6, x7, x8}

Anotherway to get primitive symbols is Unique[x]

Now we build a Groebner basis

      Timing[bs=GroebnerBasis[qs2, vars];]
       {3.14063,Null}
 Timing[sols=Solve[  (0==#&)/@bs,vars]; ]
   {3.29688,Null}

and find 2 real solutions

     relsol = Cases[strt,  
                _?(Simplify@And[(vars == Conjugate[vars]) /. #]
  qs2 /. relsol

  {{0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0}}

   {{x1->1/24,x2->1/24,x3->-(1/24),x4->1/24,x5->1/24,x6->1/24,
      x7->1/24,x8->-(1/24)},
      {x1->-(1/120),x2->-(1/120),x3->-(1/120),x4->1/24,
         x5->1/24,x6->1/24,x7->1/24,x8->-(1/24)}}

In the case of inexact coefficients either Rationalize and generate integer coefficient polynomials. Or the more general way, that is insensitive wrt. coding sloppyness, replace all numeric equality checks by tests of smallness of the abolute differences.