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 :

$$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} 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_{2} + c_{4} {\left(x_{2} + x_{3}\right)} \sqrt{\frac{{\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}}{c_{15} x_{3} + c_{0}}}}{x_{2} + x_{3}} = 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$$

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 :

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()).