Partial answer :

We firs eliminate radicals, by moving terms and squaring expressions where they are isolated. This is a "manual" task, for which replacing numerical constants by symbolic ones is of great help.

A helper function, collecting the numerical constants of a symbolic expression :

def NumConstants(Expr):
    if Expr.operator() is None:
        if Expr.is_numeric():
            return {Expr}
        return {}
    return set(reduce(union, map(NumConstants, Expr.operands())))

We use it to establish a numerical->symbolic dictionary:

NC=set(reduce(union, map(NumConstants, [eq1, eq2, eq3]))) - {1/2, -1}
SC=[var("c_{}".format(u)) for u in (1..len(NC))]
Dic=dict(zip(list(NC), SC))

As well as the reverse dictionary (using rational representation of our constants) :

Cid={Dic[u]:QQ(u) for u in Dic.keys()}

Apply this to our equations ; for this task, both Maxima's specialized methods and Sympy's simplify() are of great help :

Sys=[u.subs(Dic)==0 for u in [eq1, eq2, eq3]]
E1=(Sys[0]-Sys[0].lhs().operands()[1])^2
import sympy
E2=sympy.simplify((Sys[1] -Sys[1].maxima_methods().part(1,1,2)/Sys[1].maxima_methods().part(1,2))^2).simplify()._sage_()
E3=sympy.sympify((Sys[2]-(Sys[2].maxima_methods().part(1,1,1)/Sys[2].maxima_methods().part(1,2)))^2).simplify()._sage_()

We have now a system of symbolic equations, whose solutions contain the solution of our original system : $$\frac{{\left({\left(c_{6} x_{1} + c_{8} x_{2} + c_{8} x_{3} + c_{11}\right)} {\left(c_{15} x_{3} + c_{1}\right)} + c_{17} x_{2}\right)} c_{14}^{2} x_{1}}{c_{15} x_{3} + c_{1}} = c_{2}^{2}$$

$$\frac{{\left({\left(c_{6} x_{1} + c_{8} x_{2} + c_{8} x_{3} + c_{11}\right)} {\left(c_{15} x_{3} + c_{1}\right)} + c_{17} x_{2}\right)} c_{4}^{2} x_{1}}{c_{15} x_{3} + c_{1}} = \frac{{\left(c_{17} x_{2} + c_{9} x_{3} + c_{13}\right)} {\left(c_{5} x_{2} + c_{5} x_{3} + c_{3}\right)} c_{7}^{2} x_{2}^{2}}{{\left(c_{15} x_{3} + c_{1}\right)} {\left(x_{2} + x_{3}\right)}^{2}}$$

$$\frac{{\left({\left(c_{16} x_{2} + c_{16} x_{3} + c_{12}\right)} {\left(c_{15} x_{3} + c_{1}\right)} + c_{10} x_{2}\right)}^{2}}{{\left(c_{15} x_{3} + c_{1}\right)}^{2}} = \frac{{\left(c_{17} x_{2} + c_{9} x_{3} + c_{13}\right)} {\left(c_{5} x_{2} + c_{5} x_{3} + c_{3}\right)} c_{7}^{2} x_{3}^{2}}{{\left(c_{15} x_{3} + c_{1}\right)} {\left(x_{2} + x_{3}\right)}^{2}}$$

We can now convert these equations to a system of rational expressions, replace the symbolic constants by their values and search the roots of the numerators, keeping in mind to eliminate the denominators's roots :

R1=PolynomialRing(QQ,"x1, x2, x3")
SP=[R1((u.lhs()-u.rhs()).subs(Cid).numerator()) for u in [E1,E2,E3]]
J=R1.ideal(*SP)

But here's a snag :

sage: J.dimension()
1

Ouch ! an infinity of solutions. Two possibilities :

• our transformations may have inadvertently introduced an indeterminacy ;
• our original system has a bona fide (simple) infinity of solutions.

Of note :

sage: [{v:u.degree(v) for v in u.variables()} for u in SP]
[{x1: 2, x2: 1, x3: 2}, {x1: 2, x2: 4, x3: 4}, {x2: 4, x3: 6}]

In any case, we need something more refined than J.variety(QQbar) to get our solution...

Partial answer :

We firs eliminate radicals, by moving terms and squaring expressions where they are isolated. This is a "manual" task, for which replacing numerical constants by symbolic ones is of great help.

A helper function, collecting the numerical constants of a symbolic expression :

def NumConstants(Expr):
    if Expr.operator() is None:
        if Expr.is_numeric():
            return {Expr}
        return {}
    return set(reduce(union, map(NumConstants, Expr.operands())))

We use it to establish a numerical->symbolic dictionary:

NC=set(reduce(union, map(NumConstants, [eq1, eq2, eq3]))) - {1/2, -1}
SC=[var("c_{}".format(u)) for u in (1..len(NC))]
Dic=dict(zip(list(NC), SC))

As well as the reverse dictionary (using rational representation of our constants) :

Cid={Dic[u]:QQ(u) for u in Dic.keys()}

Apply this to our equations ; for this task, both Maxima's specialized methods and Sympy's simplify() are of great help :

Sys=[u.subs(Dic)==0 for u in [eq1, eq2, eq3]]
E1=(Sys[0]-Sys[0].lhs().operands()[1])^2
import sympy
E2=sympy.simplify((Sys[1] -Sys[1].maxima_methods().part(1,1,2)/Sys[1].maxima_methods().part(1,2))^2).simplify()._sage_()
E3=sympy.sympify((Sys[2]-(Sys[2].maxima_methods().part(1,1,1)/Sys[2].maxima_methods().part(1,2)))^2).simplify()._sage_()

We have now a system of symbolic equations, whose solutions contain the solution of our original system : $$\frac{{\left({\left(c_{6} x_{1} + c_{8} x_{2} + c_{8} x_{3} + c_{11}\right)} {\left(c_{15} x_{3} + c_{1}\right)} + c_{17} x_{2}\right)} c_{14}^{2} x_{1}}{c_{15} x_{3} + c_{1}} = c_{2}^{2}$$

$$\frac{{\left({\left(c_{6} x_{1} + c_{8} x_{2} + c_{8} x_{3} + c_{11}\right)} {\left(c_{15} x_{3} + c_{1}\right)} + c_{17} x_{2}\right)} c_{4}^{2} x_{1}}{c_{15} x_{3} + c_{1}} = \frac{{\left(c_{17} x_{2} + c_{9} x_{3} + c_{13}\right)} {\left(c_{5} x_{2} + c_{5} x_{3} + c_{3}\right)} c_{7}^{2} x_{2}^{2}}{{\left(c_{15} x_{3} + c_{1}\right)} {\left(x_{2} + x_{3}\right)}^{2}}$$

$$\frac{{\left({\left(c_{16} x_{2} + c_{16} x_{3} + c_{12}\right)} {\left(c_{15} x_{3} + c_{1}\right)} + c_{10} x_{2}\right)}^{2}}{{\left(c_{15} x_{3} + c_{1}\right)}^{2}} = \frac{{\left(c_{17} x_{2} + c_{9} x_{3} + c_{13}\right)} {\left(c_{5} x_{2} + c_{5} x_{3} + c_{3}\right)} c_{7}^{2} x_{3}^{2}}{{\left(c_{15} x_{3} + c_{1}\right)} {\left(x_{2} + x_{3}\right)}^{2}}$$

We can now convert these equations to a system of rational expressions, replace the symbolic constants by their values and search the roots of the numerators, keeping in mind to eliminate the denominators's roots :

R1=PolynomialRing(QQ,"x1, x2, x3")
SP=[R1((u.lhs()-u.rhs()).subs(Cid).numerator()) for u in [E1,E2,E3]]
J=R1.ideal(*SP)

But here's a snag :

sage: J.dimension()
1

Ouch ! an infinity of solutions. Two possibilities :

• our transformations may have inadvertently introduced an indeterminacy ;
• our original system has a bona fide (simple) infinity of solutions.

Of note :

sage: [{v:u.degree(v) for v in u.variables()} for u in SP]
[{x1: 2, x2: 1, x3: 2}, {x1: 2, x2: 4, x3: 4}, {x2: 4, x3: 6}]

In any case, we need something more refined than J.variety(QQbar) to get our solution...

sage: SP2=[R1((u.lhs().numerator()*u.rhs().denominator()-u.rhs().numerator()*u.lhs().denominator()).subs(Cid)) for u in [E1,E2,E3]]
sage: J2=R1.ideal(*SP2)
sage: J2.dimension()
1
sage: [{v:u.degree(v) for v in u.variables()} for u in SP2]
[{x1: 2, x2: 1, x3: 2}, {x1: 2, x2: 4, x3: 5}, {x2: 4, x3: 7}]

Partial answer :

We firs eliminate radicals, by moving terms and squaring expressions where they are isolated. This is a "manual" task, for which replacing numerical constants by symbolic ones is of great help.

A helper function, collecting the numerical constants of a symbolic expression :

def NumConstants(Expr):
    if Expr.operator() is None:
        if Expr.is_numeric():
            return {Expr}
        return {}
    return set(reduce(union, map(NumConstants, Expr.operands())))

We use it to establish a numerical->symbolic dictionary:

NC=set(reduce(union, map(NumConstants, [eq1, eq2, eq3]))) - {1/2, -1}
SC=[var("c_{}".format(u)) for u in (1..len(NC))]
Dic=dict(zip(list(NC), SC))

As well as the reverse dictionary (using rational representation of our constants) :

Cid={Dic[u]:QQ(u) for u in Dic.keys()}

Apply this to our equations ; for this task, both Maxima's specialized methods and Sympy's simplify() are of great help help, as well as visualizing intermediate results in "human-readable" (i. e. LaTeX) form :

Sys=[u.subs(Dic)==0 for u in [eq1, eq2, eq3]]
E1=(Sys[0]-Sys[0].lhs().operands()[1])^2
import sympy
E2=sympy.simplify((Sys[1] -Sys[1].maxima_methods().part(1,1,2)/Sys[1].maxima_methods().part(1,2))^2).simplify()._sage_()
E3=sympy.sympify((Sys[2]-(Sys[2].maxima_methods().part(1,1,1)/Sys[2].maxima_methods().part(1,2)))^2).simplify()._sage_()

We have now a system of symbolic equations, whose solutions contain the solution of our original system : $$\frac{{\left({\left(c_{6} x_{1} + c_{8} x_{2} + c_{8} x_{3} + c_{11}\right)} {\left(c_{15} x_{3} + c_{1}\right)} + c_{17} x_{2}\right)} c_{14}^{2} x_{1}}{c_{15} x_{3} + c_{1}} = c_{2}^{2}$$

$$\frac{{\left({\left(c_{6} x_{1} + c_{8} x_{2} + c_{8} x_{3} + c_{11}\right)} {\left(c_{15} x_{3} + c_{1}\right)} + c_{17} x_{2}\right)} c_{4}^{2} x_{1}}{c_{15} x_{3} + c_{1}} = \frac{{\left(c_{17} x_{2} + c_{9} x_{3} + c_{13}\right)} {\left(c_{5} x_{2} + c_{5} x_{3} + c_{3}\right)} c_{7}^{2} x_{2}^{2}}{{\left(c_{15} x_{3} + c_{1}\right)} {\left(x_{2} + x_{3}\right)}^{2}}$$

$$\frac{{\left({\left(c_{16} x_{2} + c_{16} x_{3} + c_{12}\right)} {\left(c_{15} x_{3} + c_{1}\right)} + c_{10} x_{2}\right)}^{2}}{{\left(c_{15} x_{3} + c_{1}\right)}^{2}} = \frac{{\left(c_{17} x_{2} + c_{9} x_{3} + c_{13}\right)} {\left(c_{5} x_{2} + c_{5} x_{3} + c_{3}\right)} c_{7}^{2} x_{3}^{2}}{{\left(c_{15} x_{3} + c_{1}\right)} {\left(x_{2} + x_{3}\right)}^{2}}$$

We can now convert these equations to a system of rational expressions, replace the symbolic constants by their values and search the roots of the numerators, keeping in mind to eliminate the denominators's roots :

R1=PolynomialRing(QQ,"x1, x2, x3")
SP=[R1((u.lhs()-u.rhs()).subs(Cid).numerator()) for u in [E1,E2,E3]]
J=R1.ideal(*SP)

But here's a snag :

sage: J.dimension()
1

Ouch ! an infinity of solutions. Two possibilities :

• our transformations may have inadvertently introduced an indeterminacy ;
• our original system has a bona fide (simple) infinity of solutions.

Of note :

sage: [{v:u.degree(v) for v in u.variables()} for u in SP]
[{x1: 2, x2: 1, x3: 2}, {x1: 2, x2: 4, x3: 4}, {x2: 4, x3: 6}]

In any case, we need something more refined than J.variety(QQbar) to get our solution...

EDIT :

sage: SP2=[R1((u.lhs().numerator()*u.rhs().denominator()-u.rhs().numerator()*u.lhs().denominator()).subs(Cid)) for u in [E1,E2,E3]]
sage: J2=R1.ideal(*SP2)
sage: J2.dimension()
1
sage: [{v:u.degree(v) for v in u.variables()} for u in SP2]
[{x1: 2, x2: 1, x3: 2}, {x1: 2, x2: 4, x3: 5}, {x2: 4, x3: 7}]

Same difference...

Partial answer :

We firs eliminate radicals, by moving terms and squaring expressions where they are isolated. This is a "manual" task, for which replacing numerical constants by symbolic ones is of great help.

A helper function, collecting the numerical constants of a symbolic expression :

def NumConstants(Expr):
    if Expr.operator() is None:
        if Expr.is_numeric():
            return {Expr}
        return {}
    return set(reduce(union, map(NumConstants, Expr.operands())))

We use it to establish a numerical->symbolic dictionary:

NC=set(reduce(union, map(NumConstants, [eq1, eq2, eq3]))) - {1/2, -1}
SC=[var("c_{}".format(u)) for u in (1..len(NC))]
Dic=dict(zip(list(NC), SC))

As well as the reverse dictionary (using rational representation of our constants) :

Cid={Dic[u]:QQ(u) for u in Dic.keys()}

Apply this to our equations ; for this task, both Maxima's specialized methods and Sympy's simplify() are of great help, as well as visualizing intermediate results in "human-readable" (i. e. LaTeX) form :

Sys=[u.subs(Dic)==0 for u in [eq1, eq2, eq3]]
E1=(Sys[0]-Sys[0].lhs().operands()[1])^2
import sympy
E2=sympy.simplify((Sys[1] -Sys[1].maxima_methods().part(1,1,2)/Sys[1].maxima_methods().part(1,2))^2).simplify()._sage_()
E3=sympy.sympify((Sys[2]-(Sys[2].maxima_methods().part(1,1,1)/Sys[2].maxima_methods().part(1,2)))^2).simplify()._sage_()

We have now a system of symbolic equations, whose solutions contain the solution of our original system : $$\frac{{\left({\left(c_{6} x_{1} + c_{8} x_{2} + c_{8} x_{3} + c_{11}\right)} {\left(c_{15} x_{3} + c_{1}\right)} + c_{17} x_{2}\right)} c_{14}^{2} x_{1}}{c_{15} x_{3} + c_{1}} = c_{2}^{2}$$

$$\frac{{\left({\left(c_{6} x_{1} + c_{8} x_{2} + c_{8} x_{3} + c_{11}\right)} {\left(c_{15} x_{3} + c_{1}\right)} + c_{17} x_{2}\right)} c_{4}^{2} x_{1}}{c_{15} x_{3} + c_{1}} = \frac{{\left(c_{17} x_{2} + c_{9} x_{3} + c_{13}\right)} {\left(c_{5} x_{2} + c_{5} x_{3} + c_{3}\right)} c_{7}^{2} x_{2}^{2}}{{\left(c_{15} x_{3} + c_{1}\right)} {\left(x_{2} + x_{3}\right)}^{2}}$$

$$\frac{{\left({\left(c_{16} x_{2} + c_{16} x_{3} + c_{12}\right)} {\left(c_{15} x_{3} + c_{1}\right)} + c_{10} x_{2}\right)}^{2}}{{\left(c_{15} x_{3} + c_{1}\right)}^{2}} = \frac{{\left(c_{17} x_{2} + c_{9} x_{3} + c_{13}\right)} {\left(c_{5} x_{2} + c_{5} x_{3} + c_{3}\right)} c_{7}^{2} x_{3}^{2}}{{\left(c_{15} x_{3} + c_{1}\right)} {\left(x_{2} + x_{3}\right)}^{2}}$$

We can now convert these equations to a system of rational expressions, replace the symbolic constants by their values and search the roots of the numerators, keeping in mind to eliminate the denominators's roots :

R1=PolynomialRing(QQ,"x1, x2, x3")
SP=[R1((u.lhs()-u.rhs()).subs(Cid).numerator()) for u in [E1,E2,E3]]
J=R1.ideal(*SP)

But here's a snag :

sage: J.dimension()
1

Ouch ! an infinity of solutions. Two possibilities :

• our transformations may have inadvertently introduced an indeterminacy ;
• our original system has a bona fide (simple) infinity of solutions.

Of note :

sage: [{v:u.degree(v) for v in u.variables()} for u in SP]
[{x1: 2, x2: 1, x3: 2}, {x1: 2, x2: 4, x3: 4}, {x2: 4, x3: 6}]

In any case, we need something more refined than J.variety(QQbar) to get our solution...

EDIT :

sage: SP2=[R1((u.lhs().numerator()*u.rhs().denominator()-u.rhs().numerator()*u.lhs().denominator()).subs(Cid)) for u in [E1,E2,E3]]
sage: J2=R1.ideal(*SP2)
sage: J2.dimension()
1
sage: [{v:u.degree(v) for v in u.variables()} for u in SP2]
[{x1: 2, x2: 1, x3: 2}, {x1: 2, x2: 4, x3: 5}, {x2: 4, x3: 7}]

Same difference...

PartialNOTE : answer :

We firs eliminate radicals, by moving terms and squaring expressions where they are isolated. This is a "manual" task, for which replacing totally re-written answer : the previous version were garbled by too much edits and possibly a typo-induced error.

Original data taken from the question

Check the numerical answers:

Sage

# Compute
%time SolsSM=solve([eq1, eq2, eq3], [x1, x2, x3], solution_dict=True)
print(SolsSM)
# Check
CheckSM=[[abs(v.subs(u)) for v in [eq1, eq2, eq3]] for u in SolsSM]
print(CheckSM)

CPU times: user 2min 41s, sys: 1.26 s, total: 2min 42s # Not *that* fast...
Wall time: 2min 4s
[{x1: -321.3551367159828 - 5445.311876550251*I, x2: -179.359654094372 - 223.5730259195566*I,
  x3: 2381.999033439763 - 21159.38999334523*I},
 {x1: 6.685918535164614 - 113.2918302698741*I, x2: -179.359654094372 - 223.5730259195566*I,
  x3: 2381.999033439763 - 21159.38999334523*I},
 {x1: 6.685918535164623 + 113.2918302698741*I, x2: -179.359654094372 + 223.5730259195566*I,
  x3: 2381.999033439763 + 21159.38999334523*I},
 {x1: -321.3551367159836 + 5445.311876550251*I, x2: -179.359654094372 + 223.5730259195566*I,
  x3: 2381.999033439763 + 21159.38999334523*I}]

[[2.9139133938393855e-15, 2.853183107612774e-15, 9.826835244968086e-15],
 [2.220446049250313e-16, 2.1154867529750935e-17, 9.826835244968086e-15],
 [2.344603589928194e-16, 5.288716882437734e-17, 9.826835244968086e-15],
 [3.3513599103077758e-15, 3.422400574978585e-15, 9.826835244968086e-15]]

Close enough for government's work...

Matlab

Not having Matlab, I copied the question's data and checked them :

# Solution Matlab, telle qu'imprimée dans l'énoncé
x1bar=[ -0.0116 + 0.0000*I,  0.5338 + 0.0000*I, -0.0232 + 0.0000*I,
        0.2669 + 0.0000*I, -0.0047 + 0.0000*I,  1.3191 + 0.0000*I,
        -0.0321 - 0.5445*I,  0.0007 - 0.0113*I, -0.0321 + 0.5445*I,
        0.0007 + 0.0113*I, -0.0012 + 0.0000*I,  5.1569 + 0.0000*I]
x2bar=[  1.9653 + 0.0000*I,  1.9653 + 0.0000*I,  0.2480 + 0.0000*I,
         0.2480 + 0.0000*I, -0.6495 + 0.0000*I, -0.6495 + 0.0000*I,
         -0.1794 + 0.2236*I, -0.1794 + 0.2236*I, -0.1794 - 0.2236*I,
         -0.1794 - 0.2236*I, -0.0393 + 0.0000*I, -0.0393 + 0.0000*I]
x3bar=[ -0.0174 + 0.0000*I, -0.0174 + 0.0000*I,  0.6522 + 0.0000*I,
        0.6522 + 0.0000*I,  4.8817 + 0.0000*I,  4.8817 + 0.0000*I,
        0.2382 + 2.1159*I,  0.2382 + 2.1159*I,  0.2382 - 2.1159*I,
        0.2382 - 2.1159*I,  1.7273 + 0.0000*I,  1.7273 + 0.0000*I]
SolsML=map(lambda u,v,w:{x1:1.0e4*u, x2:1.0e3*v, x3:1.0e4*w},
         x1bar, x2bar, x3bar)
# Numerical check
CheckML=[[abs(v.subs(u)) for v in [eq1, eq2, eq3]] for u in SolsML]
print(CheckML)

[[0.000102245046124216, 0.0504273057558141, 0.00407651153817091],
[0.000145426550936412, 0.0501796341587534, 0.00407651153817091],
[0.000167832577447902, 0.0000380920578297709, 0.00350920387509823]
[0.000718604724398997, 0.000924529359676793, 0.00350920387509823],
[0.000590312714312535, 1.60057187570945, 120.256029930173],
[0.00152913372345476, 1.60151069671859, 120.256029930173],
[7.89390598709807, 7.89386741627798, 0.00223166033868648],
[1.10797685571357, 1.10813035022203, 0.00223166033868648],
[7.89390598709807, 7.89386741627798, 0.00223166033868648],
[1.10797685571357, 1.10813035022203, 0.00223166033868648],
[0.000561991548113761, 1.59838923451504, 702.100941033462],
[1.06578983072011, 2.66474105678326, 702.100941033462]]

The first four answers might be correct (the original solutions were printed with a limited precision) ; the eight last seem (widely) off the mark.

Look for an "exact" answer

"Exact" has a special meaning here : the original question involves (a lot of) numerical floating-point constants, presumably approximations of some (unspecified) quantities. An "exact" solution implies treating them as "exact" decimal numbers, i. e.converting them to rationals.

Transform the equations in rational fractions

Improve readibility

In order to use Sage's polynomial handling, we have to convert those expressions to rational fractions, by isolating the terms involving square roots and squaring the equations. This is not an automated process ; in order to understand the structure of the expression, it is convenient to replace (temporarily) the numerical constants by symbolic ones is symbols :

# Utility : get the set of great help.
 
A helper function, collecting the numerical constants occurring in an expression
def GetNums(expr):
    if expr.operator() is None:
        if expr.is_numeric(): return set([expr])
        return set()
    return reduce(union, list(map(GetNums, expr.operands())))
# Set of the numerivcal constants in the three original equations
Consts=reduce(union, map(lambda u:set(GetNums(u)), [eq1, eq2, eq3]))
# Generate a set of symbolic expression constants...
SymConsts=[var("c_{}".format(u)) for u in range(len(Consts))]
# ... and build a dictionary
Dic=dict(zip(Consts,SymConsts))
# We do not want to loose the signs (represented by multiplication by -1)...
Dic.pop(-1)
# ...nor the square roots (represented by exponentiation to 1/2)
Dic.pop(1/2)

Now we can represent our original equations (without right-and side) as symbolic equations, much more readable :

def NumConstants(Expr):
    if Expr.operator() is None:
        if Expr.is_numeric():
            return {Expr}
        return {}
    return set(reduce(union, map(NumConstants, Expr.operands())))
Sys=[u.subs(Dic)==0 for u in [eq1, eq2, eq3]]

We use it to establish a numerical->symbolic dictionary:

NC=set(reduce(union, map(NumConstants, [eq1, eq2, eq3]))) - {1/2, -1}
SC=[var("c_{}".format(u)) for u in (1..len(NC))]
Dic=dict(zip(list(NC), SC))

As well as the reverse dictionary (using rational representation of our constants) which gives us the following system :

Cid={Dic[u]:QQ(u) for u in Dic.keys()}

Sys=[u.subs(Dic)==0 for u in [eq1, eq2, eq3]]
E1=(Sys[0]-Sys[0].lhs().operands()[1])^2
import sympy
E2=sympy.simplify((Sys[1] -Sys[1].maxima_methods().part(1,1,2)/Sys[1].maxima_methods().part(1,2))^2).simplify()._sage_()
E3=sympy.sympify((Sys[2]-(Sys[2].maxima_methods().part(1,1,1)/Sys[2].maxima_methods().part(1,2)))^2).simplify()._sage_()
E1=(Sys[0]-Sys[0].maxima_methods().part(1,2))^2

E2=(sympy.sympify(Sys[1]-Sys[1].maxima_methods().part(1,1,1)/\
                   Sys[1].maxima_methods().part(1,2)).simplify()._sage_())^2

E3=(sympy.sympify(Sys[2]-Sys[2].maxima_methods().part(1,1,1)/\
                    Sys[2].maxima_methods().part(1,2)).simplify()._sage_())^2

Sp=[sympy.sympify((u.lhs()-u.rhs()).factor()).simplify()._sage_().numerator() for u in [E1, E2, E3]]

R1=PolynomialRing(QQ,"x1, # Solution ring
R1=PolynomialRing(QQbar,"x1, x2, x3")
SP=[R1((u.lhs()-u.rhs()).subs(Cid).numerator()) # Translation of our symbolic constants to *rationals*
Qid={Dic[u]:QQbar(QQ(u)) for u in [E1,E2,E3]]
Dic.keys()}
SP=[R1(u.subs(Qid)) for u in Sp]
# Ideal generated by our system.
J=R1.ideal(*SP)

J.dimension()
0

sage: J.dimension()
1
J.variety(QQbar)

• our Check for null denominators for all candidate solutions (see above...).
• Check the proposed solutions against the original equations.

Our transformations may have inadvertently introduced an indeterminacy ;

Of note :

these equations where the numerical constants are replaced by their rational transcriptions, thus (theoretically) allowing to use sage: [{v:u.degree(v) for v in u.variables()} for u in SP]
[{x1: 2, x2: 1, x3: 2}, {x1: 2, x2: 4, x3: 4}, {x2: 4, x3: 6}]
is_null()

In any case, we need something more refined than J.variety(QQbar) to get our solution...

EDIT :

sage: SP2=[R1((u.lhs().numerator()*u.rhs().denominator()-u.rhs().numerator()*u.lhs().denominator()).subs(Cid)) for u in [E1,E2,E3]]
sage: J2=R1.ideal(*SP2)
sage: J2.dimension()
1
sage: [{v:u.degree(v) for v in u.variables()} for u in SP2]
[{x1: 2, x2: 1, x3: 2}, {x1: 2, x2: 4, x3: 5}, {x2: 4, x3: 7}]

Same difference...

EDIT :

Eliminating x1 between E1 and E2, and solving the result and E3 for x2 and x3 doesn't work, either...

NOTE : This is a totally re-written answer : the previous version were garbled by too much edits and possibly a typo-induced error.

Original data taken from the question

Check the numerical answers:

Sage

# Compute
%time SolsSM=solve([eq1, eq2, eq3], [x1, x2, x3], solution_dict=True)
print(SolsSM)
# Check
CheckSM=[[abs(v.subs(u)) for v in [eq1, eq2, eq3]] for u in SolsSM]
print(CheckSM)

CPU times: user 2min 41s, sys: 1.26 s, total: 2min 42s # Not *that* fast...
Wall time: 2min 4s
[{x1: -321.3551367159828 - 5445.311876550251*I, x2: -179.359654094372 - 223.5730259195566*I,
  x3: 2381.999033439763 - 21159.38999334523*I},
 {x1: 6.685918535164614 - 113.2918302698741*I, x2: -179.359654094372 - 223.5730259195566*I,
  x3: 2381.999033439763 - 21159.38999334523*I},
 {x1: 6.685918535164623 + 113.2918302698741*I, x2: -179.359654094372 + 223.5730259195566*I,
  x3: 2381.999033439763 + 21159.38999334523*I},
 {x1: -321.3551367159836 + 5445.311876550251*I, x2: -179.359654094372 + 223.5730259195566*I,
  x3: 2381.999033439763 + 21159.38999334523*I}]

[[2.9139133938393855e-15, 2.853183107612774e-15, 9.826835244968086e-15],
 [2.220446049250313e-16, 2.1154867529750935e-17, 9.826835244968086e-15],
 [2.344603589928194e-16, 5.288716882437734e-17, 9.826835244968086e-15],
 [3.3513599103077758e-15, 3.422400574978585e-15, 9.826835244968086e-15]]

Close enough for government's work...

Matlab

Not having Matlab, I copied the question's data and checked them :

# Solution Matlab, telle qu'imprimée dans l'énoncé
x1bar=[ -0.0116 + 0.0000*I,  0.5338 + 0.0000*I, -0.0232 + 0.0000*I,
        0.2669 + 0.0000*I, -0.0047 + 0.0000*I,  1.3191 + 0.0000*I,
        -0.0321 - 0.5445*I,  0.0007 - 0.0113*I, -0.0321 + 0.5445*I,
        0.0007 + 0.0113*I, -0.0012 + 0.0000*I,  5.1569 + 0.0000*I]
x2bar=[  1.9653 + 0.0000*I,  1.9653 + 0.0000*I,  0.2480 + 0.0000*I,
         0.2480 + 0.0000*I, -0.6495 + 0.0000*I, -0.6495 + 0.0000*I,
         -0.1794 + 0.2236*I, -0.1794 + 0.2236*I, -0.1794 - 0.2236*I,
         -0.1794 - 0.2236*I, -0.0393 + 0.0000*I, -0.0393 + 0.0000*I]
x3bar=[ -0.0174 + 0.0000*I, -0.0174 + 0.0000*I,  0.6522 + 0.0000*I,
        0.6522 + 0.0000*I,  4.8817 + 0.0000*I,  4.8817 + 0.0000*I,
        0.2382 + 2.1159*I,  0.2382 + 2.1159*I,  0.2382 - 2.1159*I,
        0.2382 - 2.1159*I,  1.7273 + 0.0000*I,  1.7273 + 0.0000*I]
SolsML=map(lambda u,v,w:{x1:1.0e4*u, x2:1.0e3*v, x3:1.0e4*w},
         x1bar, x2bar, x3bar)
# Numerical check
CheckML=[[abs(v.subs(u)) for v in [eq1, eq2, eq3]] for u in SolsML]
print(CheckML)

[[0.000102245046124216, 0.0504273057558141, 0.00407651153817091],
[0.000145426550936412, 0.0501796341587534, 0.00407651153817091],
[0.000167832577447902, 0.0000380920578297709, 0.00350920387509823]
[0.000718604724398997, 0.000924529359676793, 0.00350920387509823],
[0.000590312714312535, 1.60057187570945, 120.256029930173],
[0.00152913372345476, 1.60151069671859, 120.256029930173],
[7.89390598709807, 7.89386741627798, 0.00223166033868648],
[1.10797685571357, 1.10813035022203, 0.00223166033868648],
[7.89390598709807, 7.89386741627798, 0.00223166033868648],
[1.10797685571357, 1.10813035022203, 0.00223166033868648],
[0.000561991548113761, 1.59838923451504, 702.100941033462],
[1.06578983072011, 2.66474105678326, 702.100941033462]]

The first four answers might be correct (the original solutions were printed with a limited precision) ; the eight last seem (widely) off the mark.

Look for an "exact" answer

"Exact" has a special meaning here : the original question involves (a lot of) numerical floating-point constants, presumably approximations of some (unspecified) quantities. An "exact" solution implies treating them as "exact" decimal numbers, i. e.converting them to rationals.

Transform the equations in rational fractions

Improve readibility

In order to use Sage's polynomial handling, we have to convert those expressions to rational fractions, by isolating the terms involving square roots and squaring the equations. This is not an automated process ; in order to understand the structure of the expression, it is convenient to replace (temporarily) the numerical constants by symbols :

# Utility : get the set of numerical constants occurring in an expression
def GetNums(expr):
    if expr.operator() is None:
        if expr.is_numeric(): return set([expr])
        return set()
    return reduce(union, list(map(GetNums, expr.operands())))
# Set of the numerivcal constants in the three original equations
Consts=reduce(union, map(lambda u:set(GetNums(u)), [eq1, eq2, eq3]))
# Generate a set of symbolic constants...
SymConsts=[var("c_{}".format(u)) for u in range(len(Consts))]
# ... and build a dictionary
Dic=dict(zip(Consts,SymConsts))
# We do not want to loose the signs (represented by multiplication by -1)...
Dic.pop(-1)
# ...nor the square roots (represented by exponentiation to 1/2)
Dic.pop(1/2)

Now we can represent our original equations (without right-and side) as symbolic equations, much more readable :

Sys=[u.subs(Dic)==0 for u in [eq1, eq2, eq3]]

which gives us the following system :

$$\begin{align} c_{14} $$c_{14} \sqrt{\frac{{\left({\left(c_{16} x_{1} + c_{7} x_{2} + c_{7} x_{3} + c_{10}\right)} {\left(c_{15} x_{3} + c_{0}\right)} + c_{18} x_{2}\right)} x_{1}}{c_{15} x_{3} + c_{0}}} + c_{1} & = 0 \\ \frac{\sqrt{\frac{{\left(c_{18} 0$$

$$\frac{{\left(c_{15} x_{3} + c_{0}\right)} \sqrt{\frac{{\left(c_{18} x_{2} + c_{8} x_{3} + c_{13}\right)} {\left(c_{5} x_{2} + c_{5} x_{3} + c_{3}\right)}}{c_{15} x_{3} + c_{0}}} c_{6} x_{3} + {\left(c_{12} x_{2} + c_{12} x_{3} + c_{11}\right)} {\left(c_{15} x_{3} + c_{0}\right)} {\left(x_{2} + x_{3}\right)} + c_{9} x_{2} + c_{9} x_{3}}{{\left(c_{15} x_{3} + c_{0}\right)} {\left(x_{2} + x_{3}\right)}} & = 0 \end{align}$$0$$

Eliminate radicals

For this task, two "external" helpers are really useful :

• Maxima's part function is a great shortcut for expression surgery ;
• Sympy's simplify is often the fastest way to a concise rewriting of a complicated expression.

So :

import sympy
E1=(Sys[0]-Sys[0].maxima_methods().part(1,2))^2

$$\frac{{\left({\left(c_{16} x_{1} + c_{7} x_{2} + c_{7} x_{3} + c_{10}\right)} {\left(c_{15} x_{3} + c_{0}\right)} + c_{18} x_{2}\right)} c_{14}^{2} x_{1}}{c_{15} x_{3} + c_{0}} = c_{1}^{2}$$

E2=(sympy.sympify(Sys[1]-Sys[1].maxima_methods().part(1,1,1)/\
                   Sys[1].maxima_methods().part(1,2)).simplify()._sage_())^2

$$\frac{{\left({\left(c_{16} x_{1} + c_{7} x_{2} + c_{7} x_{3} + c_{10}\right)} {\left(c_{15} x_{3} + c_{0}\right)} + c_{18} x_{2}\right)} c_{4}^{2}}{c_{15} x_{3} + c_{0}} = \frac{{\left(c_{18} x_{2} + c_{8} x_{3} + c_{13}\right)} {\left(c_{5} x_{2} + c_{5} x_{3} + c_{3}\right)} c_{6}^{2} x_{2}^{2}}{{\left(c_{15} x_{3} + c_{0}\right)} {\left(x_{2} + x_{3}\right)}^{2}}$$

E3=(sympy.sympify(Sys[2]-Sys[2].maxima_methods().part(1,1,1)/\
                    Sys[2].maxima_methods().part(1,2)).simplify()._sage_())^2

$$\frac{{\left(c_{18} x_{2} + c_{8} x_{3} + c_{13}\right)} {\left(c_{5} x_{2} + c_{5} x_{3} + c_{3}\right)} c_{6}^{2} x_{3}^{2}}{{\left(c_{15} x_{3} + c_{0}\right)} {\left(x_{2} + x_{3}\right)}^{2}} = \frac{{\left({\left(c_{12} x_{2} + c_{12} x_{3} + c_{11}\right)} {\left(c_{15} x_{3} + c_{0}\right)} {\left(x_{2} + x_{3}\right)} + c_{9} x_{2} + c_{9} x_{3}\right)}^{2}}{{\left(c_{15} x_{3} + c_{0}\right)}^{2} {\left(x_{2} + x_{3}\right)}^{2}}$$

which can be rewritten as a system of expressions (equations with implicit null right-hand-side)

Sp=[sympy.sympify((u.lhs()-u.rhs()).factor()).simplify()._sage_().numerator() for u in [E1, E2, E3]]

(The LaTeX representation is too large for the available engine...)

Solve the system of fractions

As far as I know, Sage does not have (built-in) facilities for solving rational fractions. So we will solve the fraction denominators, while keeping in mind to eliminate "solutions" where any of the denominators is null.

# Solution ring
R1=PolynomialRing(QQbar,"x1, x2, x3")
# Translation of our symbolic constants to *rationals*
Qid={Dic[u]:QQbar(QQ(u)) for u in Dic.keys()}
SP=[R1(u.subs(Qid)) for u in Sp]
# Ideal generated by our system.
J=R1.ideal(*SP)

J.dimension()
0

So far, so good... But there is a fly in the ointment :

J.variety(QQbar)

hasn't returned for > 8 hours... CPU usage shows that one process (Sage...) eats 100% of its CPU, RAM usage is absolutely stable. I begin to suspect a bug. Advice welcome...

What should follow

• Check for null denominators for all candidate solutions (see above...).
• Check the proposed solutions against the original equations.equations :

Our transformations may have introduced extraneous solutions ($x^2=a^2$ has more solutions than $x=a$). The simplest way is to check the solutions against the original equations (or, better, against a version of these equations where the numerical constants are replaced by their rational transcriptions, thus (theoretically) allowing to use is_null()).