You may have made the mistake of interpreting **x=y^(1/3)** (*outside* of the **RealDomain** context) as being purely real-valued for **y::real**.

Consider which of the results in **s1** below will evaluate to **x=-2** when **y=-8** . Notice that it is *not* the result of **x=y^(1/3)** .

restart;
s1 := [solve( {y-x^3}, x)];
[ / (1/3)\ / 1 (1/3) 1 (1/2) (1/3)\
[{ x = y }, { x = - - y + - I 3 y },
[ \ / \ 2 2 /
/ 1 (1/3) 1 (1/2) (1/3)\ ]
{ x = - - y - - I 3 y }]
\ 2 2 / ]
for s in s1 do
print( s, radnormal(eval(s, y=-8)), evalf[5](evalf(eval(s, y=-8))) );
end do:
/ (1/3)\ / (1/2)\
{ x = y }, { x = 1 + (-3) }, {x = 1.0000 + 1.7321 I}
\ / \ /
/ 1 (1/3) 1 (1/2) (1/3)\
{ x = - - y + - I 3 y }, {x = -2}, {x = -2.0001 + 0. I}
\ 2 2 /
/ 1 (1/3) 1 (1/2) (1/3)\ / (1/2)\
{ x = - - y - - I 3 y }, { x = 1 - I 3 }, {x = 1.0001 - 1.7321 I}
\ 2 2 / \ /

But you do want **x=y^(1/3)** when **y>0**. So (outside of using **RealDomain**) a purely real-valued solution might be expressed as a **piecewise**, using the first item in **s1** when **y>=0** and the second item in **s1** when **y<0** . For example,

s2 := piecewise( y>=0, y^(1/3), y<0, -(-y)^(1/3) ):
plot( s2, y=-8..8, size=[200,200] );

Does anyone see a nice way to get **:-solve** to return something like that **piecewise**?

An alternative is to get **solve** to return an implicit **RootOf** (Maple 2015.2 used here) which evaluates to a real for real numeric **y**. For example,

s3 := solve( y-x^3, x, parametric, real );
[[ / 3 \]]
[[x = RootOf\_Z - y, index = real[1]/]]
s := s3[1];
[ / 3 \]
[x = RootOf\_Z - y, index = real[1]/]
evalf(eval(s,y=-8));
[x = -2.]

Unfortunately this is slower to plot, as it evaluates to floats internally using **`evalf/RootOf`** which I suspect has to go through **fsolve** or some numeric root-finder.

plot( eval(x,s), y=-8..8, size=[200,200] );

pr.mw