I want to study an algebraic curve defined by equations of the form
$$ a_1 \sqrt{f_1(x)} + ... + a_n \sqrt{f_n(x)} = 0, $$
where $x$ is a real variable and $f_i$ are polynomials. $ a_1,... a_n $ could be functions of $x$, and there could be more than one equation (that is a system), and it has also constant parameters (not just numbers) and other variables ($y$ for two dimensions), though all this may be not important for this question.
I wanted to use SageMath, but it seems it can't even solve an equation
$$ \sqrt{x} + \sqrt{1-x} = 1 $$
with obvious solutions 0 and 1.
x = var('x')
rad_pol = sqrt(x) + sqrt(1-x) - 1
solve(rad_pol, x)
# Out: [sqrt(x) == -sqrt(-x + 1) + 1]
I checked Mathematica online, and it solves such an equation.
However, I can't be sure that I will be able to solve my systems of equations online (because of proprietary software restrictions or online limitations), and maybe there exists an even more powerful instrument to solve such equations. In fact, I'm more interested in a good instrument to solve a system of equations, because I could do radical simplifications manually, but in this respect I see little difference between systems, because the Gröbner basis is almost universally implemented. However, a decent computer algebra system may have more advanced methods to solve radical equations than I could do manually (like this method with resultants).
I checked that online Sage can neither solve the simple equation, so the problem is likely not my installation. Do you know good free software to do that? If not, what could be considered the most powerful software for these kinds of problems, even non-free?
UPD: thanks to Samuel Lelièvre's answer, we can solve the equation above in Sage. Unfortunately, it still can't solve an equation for two points in two dimensions. Suppose we want to find a solution for
$$|\mathbf{n}_1 + \mathbf{n}_2| = 1,$$
where
$$\mathbf{n}_i = \frac{(x_i - x_0; y_i - y_0)}{\sqrt{(x_i - x_0)^2 + (y_i - y_0)^2}},$$
or, by squaring the first equation, $1 + 2 \mathbf n_1 \mathbf n_2 = 0$,
$$ \sqrt{(x_1-x_0)^2 + (y_1-y_0)^2} \sqrt{(x_2-x_0)^2 + (y_2-y_0)^2} + 2 \left((x_1-x_0)(x_2-x_0) + (y_1-y_0)(y_2-y_0)\right) = 0$$
This is a curve on which the sum of unit vectors to two given points is equal to one. The answer is known to be a part of a circle, and it can be solved by Mathematica (both with numeric point coordinates and with symbolic ones). Again, it can't be solved by sagemath, unfortunately. So the question about a powerful enough system remains.
# variables
x0, y0 = var('x0 y0')
# parameters
x1, x2, y1, y2 = 0, 1, 0, 0 # var('x1, x2, y1, y2')
r1, r2 = var('r1, r2')
# can't solve
eq2 = 2 * ((x1-x0)*(x2-x0) + (y1-y0)*(y2-y0)) == - r1 * r2
solve([eq2, r1**2 - ((x1-x0)**2 + (y1-y0)**2), r2**2 - ((x2-x0)**2 + (y2-y0)**2)], x0, y0, to_poly_solve=True)
# Mathematica succeeds:
# solve 2 * ((x1-x0)*(x2-x0) + (y1-y0)*(y2-y0)) == - sqrt((x1-x0)**2 + (y1-y0)**2) * sqrt((x2-x0)**2 + (y2-y0)**2), x1=0, y1=0, x2=1, y2=0 for x0, y0
#
# maybe try without new variables for radicals?
eq3 = 2 * ((x1-x0)*(x2-x0) + (y1-y0)*(y2-y0)) + sqrt((x1-x0)**2 + (y1-y0)**2) * sqrt((x2-x0)**2 + (y2-y0)**2) == 0
solve(eq3, x0, y0, to_poly_solve=True)
# doesn't work again, empty set.