Aaah ! A nice one...

Trying to solve :

```
Sol = solve(Sys, SVars)
```

Long error message : elided, seemingly false,... **but another error message during the processing of the first :**

```
TypeError: unable to make sense of Maxima expression '[if((-pi/2 < parg(374-2*c498417)) and (-pi/2 < parg(c498417-187)) and (-pi/2 < parg(c498417)) and (parg(374-2*c498417) <== pi/2) and (parg(c498417-187) <== pi/2) and (parg(c498417) <== pi/2),[_SAGE_VAR_H == -4*c498417,_SAGE_VAR_K == ((-2*c498417^2)-70)/187,_
```

... much elided again.

It seems that Maxima, duly piloted by Sage, is able to find a solution, but Sage is unable to translate it. The error message seems to point to the test of the argument of a quantity (possibly checking for the sign of a radicand).

This second error may indicate that the firs one has been derailed, hence a dubious meaning...

Further attempts with Mathematica (or the *gratis* Wolfram Engine, used through emacs' `wolfram-mode`

) point to the same direction.

Your system has three equations involving square roots. We can try to eliminate those radicals, by replacing the equations by other ones whise roots are a *superset* of the original ones.

It happens that `sympy`

has a nice tool for this : from `sympy.solvers.solvers.unrad?`

:

```
Remove radicals with symbolic arguments and return (eq, cov), None,
or raise an error.
None is returned if there are no radicals to remove.
NotImplementedError is raised if there are radicals and they cannot
be removed or if the relationship between the original symbols and
the change of variable needed to rewrite the system as a polynomial
cannot be solved.
Otherwise the tuple, "(eq, cov)", is returned where:
*eq*, "cov"
*eq* is an equation without radicals (in the symbol(s) of
interest) whose solutions are a superset of the solutions to the
original expression. *eq* might be rewritten in terms of a new
variable; the relationship to the original variables is given by
"cov" which is a list containing "v" and "v**p - b" where "p" is
the power needed to clear the radical and "b" is the radical now
expressed as a polynomial in the symbols of interest. For
example, for sqrt(2 - x) the tuple would be "(c, c**2 - 2 + x)".
The solutions of *eq* will contain solutions to the original
equation (if there are any).
```

[ Details elided but capital ! ]

So, let's build a system whose roots are a superset of the original ones :

```
import sympy
USys=[u[1][0]._sage_() if u[1] is not None else u[0]
for u in list(zip(Sys,
map(lambda u:sympy.solvers.solvers.unrad(sympy.sympify(u)), Sys)))]
```

Use it to get a list of *candidate* solutions :

```
sage: %time USol=solve(USys, SVars)
CPU times: user 1.81 s, sys: 16 ms, total: 1.83 s
Wall time: 1.45 s
```

and keep the solutions that check the original system :

```
CUSol=[s for s in USol if all([bool(e.subs(s).simplify_full()) for e in Sys])]
sage: len(CUSol)
5
```

But all these solutions are not distinct :

```
sage: len(list(Set(set(u) for u in CUSol)))
4
```

FWIW :

Mathematica's `Solve`

proposes these 4 solutions, but from a direct interface ; Sage's interface gets derailed by Mathematica's *warning* message, interpreted by the interface as an *error* message.

Mathematica's `Reduce`

points at *possible* other solutions, which I have not (yet) checked.

Morality : equations system solving is not (yet) automatic. Sage and consorts may (seriously) *compute* for you, but won't *think* for you.

HTH,

A minimal command raising the same error seems to be:

I created https://trac.sagemath.org/ticket/34659 for the simplified issue.

Added

`confirmed_bug`

tag, see https://ask.sagemath.org/questions/sc...