Another four tries and I finally nailed it:

Explore(   
   plot3d(
      [[abs(sqrt(z^2+c)), theta, z], [abs(sqrt(z^2+c)), theta, -z]],
      z= `if`(c >= 0, -sqrt(8), sqrt(abs(c)))..sqrt(8), theta= -Pi..Pi,
      grid= [51, 36], 
      coords= cylindrical,
      #stabilize the axes' lengths for Explore 
      view= [-2..2, -2..2, -sqrt(8)..sqrt(8)],
      thickness= 0 #use thinness gridlines available
    ),
    #decimal points required for float Explore parameters,
    #lest they be treated like integers:
    c= -1.0..1.0 
);

Replace the entire for loop at the end with

(LI,LR):= selectremove(hastype, L2, nonreal(float));

Then LI will contain all the you solutions that have an imaginary component, and LR will contain the rest.

Note that real numbers are also considered to be complex, so searching for complex(...)  types is ineffective. Also, if you search with type, then the entire expression being compared will need to match the given type, which won't work because the items being checked are not numbers but sets of equations. That's why I use hastype instead. So, another way to search is

(LR,LI):= selectremove(type, L2, set(name = realcons));

To avoid the error, use add instead of sum. I'm assuming that RM is some coefficient matrix that's already defined? Then remove the unevaluation quotes.

If you continue to have trouble, please post a worksheet, including the values of RM.

I'm not sure what you mean by "associate". Do you mean that you want to make a new array composed of the two original arrays? Then, let's suppose that you have two one-dimensional arrays A and B of the same length, L. You can compose them into a 2xL array with

< A, B >

You can compose them into an Lx2 array with

< A^+ | B^+ >

The +-sign exponents mean "transpose".

There are many other ways that I can interpret "associate" that all have easy solutions in Maple. A type of association that's a bit esoteric but that I often use anyway is to define a function f which maps each element of A to its corresponding element of B:

f:= 'f': assign~(f~(A) =~ B):

You did HeffR:= convert(Heff, rational), but then you did Eigenvalues(Heff) rather than Eigenvalues(HeffR).

Assignments to sets and lists (whether by := or by subsop) are extremely inefficient. It is only begrudingly that Maple allows it for lists, and even that's only for small lists (upto size 100 I think). I think that that's allowed for the benefit of newbies doing homework problems (where efficiency doesn't matter much) who want to use lists as vectors or arrays. But nobody (in their right mind) wants to use a set as an array.

Assigning to set positions is even more inefficient than it is for lists because the new element needs to be put in the correct position (since sets are always stored sorted). So, even if S[k]:= x were allowed, after doing it, it wouldn't necessarily be true that S[k] = x or that nops(S) was the same as it was before. Those two facts alone make S[k]:= x too weird to consider using, regardless of efficiency concerns.

So, whenever you feel the urge to do S[k]:= x, you need to rethink that and come up with a better way to code it. There's always a better way. MaplePrimes can offer help in any situation like that.

If you change solve to eliminate, you'll get a less-trivial solution: a = (k^2-1)/3, b= 0, c= 0.

The Adams-Bashforth method is implemented in dsolve via

dsolve({...}, numeric, method= classical[adambash]);

See ?dsolve,numeric,classical

Use

int(1/(1+sqrt(x)), x= 0..1, method= ftoc);

To aid in remembering that, ftoc stands for fundamental theorem of calculus.

For the real-rooted polynomial, use the characteristic polynomial of a very sparse symmetric matrix all of whose entries are in {-1, 0, 1} of order 600. Then multiply that by any integer polynomial of degree 200 and expand. This only takes a few seconds (in modern Maple):

P:= expand(
   LinearAlgebra:-CharacteristicPolynomial(
      LinearAlgebra:-RandomMatrix(
         600, shape= symmetric, density= 1/200, generator= rand(-1..1)
      ), x
   )*
   randpoly(x, degree= 200)
);

Another way is to make the 1/2 a symbol, do the integration, then make the symbol equal 1/2. Thus:

restart:
f:= ((r-I*t)^(-s)-(r+I*t)^(-s))/(I/r):
int(f, t= 0..infinity) assuming s > 1, r > 0:
simplify(eval(%, r= 1/2));

On the second pass through the final loop (thus, when i=1), the dsolve doesn't find any solutions and it returns NULL. Thus S[1] is the NULL sequence. Thus S[1][1] has an "invalid subscript selector" because you're asking for the first member of an empty sequence.

I did manage to get a numeric solution from dsolve with no trouble.

If you want to apply the operation to column j of matrix A and put the results back into column j (thereby overwriting the original column j), you can do

A[..,j]:= evalf(convert~(A[..,j], temperature, Celsius, kelvin)):

Note that kelvin is with a lowercase k (because Kelvin means something else in Maple, completely unrelated to temperature).

In Maple, negative reals raised to fractional powers return the principal complex value. This is the best general answer for a system that generally handles complex numbers. But Maple also provides the command surd to handle the real case. So,

f:= x-> surd(4*x^2-4, 5)^4;

That means the 5th root to 4th power.

And here's another solution that I don't like as much, but I'll mention it because if I don't someone else likely will:

f:= proc(x) uses RealDomain; (4*x^2-4)^(4/5) end proc:

There are many ways to do it. I'd do

R:= unapply(index~(<P1, P2, P3, P4, P1 | Q1, Q2, Q3, Q4, Q1>, i), i):
plot([seq(R(i), i= 1..21)], color= plots:-setcolors()[1]);

That'll plot them all in the same color. If you'd like different colors, just remove the color option completely. The index command in the first line requires Maple version at least "Maple 2017". If you have an earlier version, I can very easily give you another command. Let me know.