@Johan159 fsolve works here.

 M_V_as_function_of_hight_cc_temporary_v3_dharr.mw

@Johan159 Sorry, I misunterstood what you are trying to do, and I agree you want to solve for y_. Look at the implicitplot here - there are two branches - which do you want (or both)?.

The problem is too hard with the double piecewise, since for a given y, the answer for y_ can be in any of three regions, so there are 9 possibilities overall that Maple needs to consider. I suggest you have separate equations for each y region,  then solve for y_ (perhaps still piecewise, but better not). Even then, there might not be an analytical solution from solve and you then have to use fsolve (and make a function that takes uses fsolve to output y_ given y.. I suugest as a starting point you look closely at the implicitplot for one of the regions you are interested in - choose a region for y that has y_ in only one region and try to solve that first (no piecewise). The you will know whether you can use solve. 

I would say in general solutions with piecewise are difficult to solve, but easier with fsolve.

M_V_as_function_of_hight_cc_temporary_dharr.mw

@Aixleft math A quick way is 

pcube := draw(poly(color = black, cutout = 0.98));

though cutouts don't really leave lines (try a smaller number to see).

Alternatively you could just find the edges and plot them as lines. Strangely, geom3d doesn't have a way to find edges directly; you have to do it via faces. See attached.

Download cube.mw

@Carl Love That's an interesting assumption about how is works. Since evalb doesn't want to work with sqrt(3) and other algebraic numbers, and is will, I assumed that is is purely symbolic, and evaluates is(sqrt(3)<3)by squaring both sides or similar, so floats are never required (and hence takes much longer). I like to think there is never any doubt when Maple works symbolically, and it is probably possible to always find an example where float numerics fail.

If I change the 2D math to 1D math and execute with the ! icon it works.

@Hullzie16 OK, solve,identity seems to work only up to about -5, and as you say doesn't take into account the fact that you can solve them sequentially. So your method seems OK to me; a slight improvement might be

for i to 10 do b[i]:=factor(solve(coeff(E2-Q/r,r,-i),b[i])) end do;

@GFY A variation of  @C_R 's changes. Not exactly sure what you want to do with the eigenvalues. Changes in red.

equation915_debug.mw

@sursumCorda Yes, thanks for pointing that out; I missed the distnctness aspect. You could separately do the cases of 1, 2 and 3 distinct roots, but I suppose postprocessing is required.

dec1 := RealRootClassification([f], [], [x], [p, q, r], 3, 1, R);
dec2 := RealRootClassification([f], [], [x], [p, q, r], 3, 2, R);
dec3 := RealRootClassification([f], [], [x], [p, q, r], 3, 3, R);

The [p,q,r] was to look only for cases where these are non-zero, but seems to make no difference here over []. However, when I was playing with RealTrangularize there was a difference between [p,q,r] and [], so I'm confused about this. 

discrim(f,x) is zero iff there exist multiple roots, so that may also help. Here it is just the now familiar

 -4*p^3*r + p^2*q^2 + 18*p*q*r - 4*q^3 - 27*r^2

and Descarte's rule of signs is useful in this case, though you did want it without polynomials theory.

@dharr This is very close, but has the condition p^2*q + 3*p*r - 4*q^2>0 instead of q<0 - perhaps they are equivalent but that wasn't immediately obvious to me.

Download cubic2.mw

Setting all variables to 1 shows that the odetest result cannot be zero

q := odetest(G, F4);
eval(q, indets(q, name) =~ 1);

Since many of your equations are set up by pasting output, it's possibly a cut and paste error, but that makes it nearly impossible to diagnose.

@Saha So you will need to make a loop changing the L values. Each time though the loop use dsolve and find Q, and store L^2 and Q in Vectors. Then you can use pointplot to plot the vector of Q values vs the vector of L^2 values.

Your worksheet doesn't have theta, so I'm guessing you want f'(1). The output of dsolve can be used to get this, e.g., c1(1) gives

[x = 1., f(x) = 1.00000000000000, diff(f(x), x) = 0.294118092731875] 

and so to extract the value 0.294118092731875 use eval(diff(f(x), x), c1(1)).

@Carl Love I added a 60 Digit calculation to the 30 digit one

evalf[30](eval([p, sinnx, p - sinnx], x = -1.95));
evalf[60](eval([p, sinnx, p - sinnx], x = -1.95));

 (errors (p-sinnx) ca. 10^(-13) for 30 digits and 10^(-43) for 60 digits.)

At 30 digits I was thinking there was a true error because Maple usually does much better than half the digits correct. But since the functions are identical, I increased to 60, which shows how poor the 30 Digit one was.

@Carl Love Thanks for the resolution of the identical/approx issue. Nice solution for the plot - vote up.

@FDS I get the same as you on Maple 2024.1 on Windows. Opening the Excel file in Excel (Office 365) everything looks fine. So that is a mystery; perhaps contact technical support.