You have many subscripted Energy variables so I can't really understand what you want to do, and some change of variables is mentioned but without being specific. It seems strange to differentiate and then integrate almost the same thing. It would be a lot easier to understand if you have variables instead of numbers (see my example below).

Perhaps you could just solve a differential equation. If you just want plots, something like the following implicitly dies a change of variables from t to E. Or you could use PDEtools:-dchange to get Maple to do the change of variables on the original ode. Perhaps the following does something that can be adapted to your situation.

Download DAE.mw

@mary120 You see at the end of the worksheet that the numerator has M^6 + const*M^4+const*M^2+const and the denominator has M^2*(M^6 +const*M^4+const*M^2+const)

@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.