@ajfriedlan See the difference in :
(1,2,3)[1], (1)[1];
[1,2,3][1], [1][1];

Your code uses the first version, it should be using the second. That is:
 Rootz := [solve](epsilon__k(q, J__a[l])-epsilon__k(kappa, J__a[l])+.16*e = 0);

You can always compose plots in the end, so I suggest you work on them sparately first.

@WJRAGG Simply evaluate the following list: [K, P, Q, ap]; and print the output here.

The parameters : {K, P, Q, ap} are missing from your code.

@Carl Love

g := subs(F = unapply(f,a,b,x), (a,b) -> (x-> F(a, b, x)));

@Carl Love This g does not do what was intended.

Also I have a Notification "Carl Love mentioned you in a Answer". Why does it say that?

@666 jvbasha or plots[pointplot3d]([seq(seq([bix[i],prx[j],X1[i,j]],i=1..NN),j=1..NN1)]);

@666 jvbasha try: print(eval(X1));

@Carl Love I dont think using indets will ever match subsindets code. I expect subsindets recurses using type & hastype.

Another difference:

e:= f(g(y), f(f(f(x), g(y)), f(x)), g(y)):
i:= 0:  SubsIndets(e, function, h);
i:= 0:  subsindets(e, function, h);

@ecterrab That was my mistake. I forgot type testing was done to the initial value not to the resulting value. "You find that which you seek".

@Carl Love The reason is subsindets uses strict operand ordering unlike the indets call used.
Try this example:

restart;
i := 0;
h := proc(x) global i; i := i + 1; T[i, op(0,x)](op(x)) end;
subsindets(f(f(f(f(x), g(y)), f(x)), g(y)), function, h);

@Carl Love I dont know why this selected the best answer. It gives undesirable elements in the result.

Because y is in your expression:

expr := y^2*sin(1/y) + y^(3/2) + y + x*y^7;

y is included in the result. Should y be included in the answer to the following expression:

expr := y^2*sin(1/y) + y^(3/2) + x*y^7;

I have noticed that after entering the login/password "My Maple Cloud" section momentarily appears and vanishes and I remain not logged in. I am able to login via web browser, but I would prefer to login via my Maple.

@9009134 Just like Carl Love pointed out, the name differs from its appearance. That is why I had to use:

S := remove(type, indets(B, function), trig);

to collect the functions. You can get the names & substitutions from this also:

S2 := (sort@[op]@map2)(op, 0, remove(has, S, diff));
S3 := [f__31, f__21, f__11];
lprint(`=`~(S2, S3));
subs(`=`~(S2, S3), eval(W));