@Jean-Michel  
The  ListTools:-Search command allows you to find the positions of only the first-level operands in a certain list.
Your example:

restart;
with(ListTools):
M := [[1, a], a, [[3 + a]], 2*a - 1]:
Search(a, M);
Search(a, M[1]);
Search(2*a-1, M);

@nm  Your method doesn't work in older versions of Maple (like my Maple 2018.2). However, the following method works successfully:

restart;
a, b, c, d := Matrix([[1,2],[3,4]]), Matrix([[1,2,3],[4,5,6],[7,8,9]]),
Matrix([[1,2,3,4],[5,6,7,8],[9,10,11,12],[13,14,15,16]]),
Matrix([[1,2,3,4,5],[6,7,8,9,10],[11,12,13,14,15],[16,17,18,19,20],[21,22,23,24,25]]);
A:=Matrix(2,2,{(1,1)=a,(1,2)=b,(2,1)=c,(2,2)=d});

@nm The code uses coloring using the options provided by  plots:-densityplot  command. Individual coloring of each line is not provided. All lines are colored with the same color (black by default). Implementing your wishes would require significant changes to the procedure code.

@Jean-Michel  Should be an ordered range, but may be  -sin(x) > cos(x) , for example -sin(-Pi/2) .. cos(-Pi/2) = 1 .. 0

Correction:

restart;
with(RandomTools):
l := x -> Generate(float(range = min(-sin(x),cos(x)) .. max(-sin(x),cos(x)))):
randomize():
plot([-sin(x),cos(x),'l(x)'], x = -Pi .. Pi);

@nm  Try just increasing those indents in  the_legend  a little. Or, while maintaining the same spaces, slightly reduce the size of the plot, for example, add the option  size=[350, 350]

The logic behind all these decisions is that both legends would not fit in length, being located on one line.

 @Matt C Anderson  
Your method for constructing the set of desired objects using  b := `union`(b, {a}) construction is extremely inefficient. This isn't particularly evident in your example of finding twin primes, as they are quite rare, but in other cases, the difference compared to other methods would be very significant.

The examples below compare the construction of the set of all odd primes less than 1,000,000 using your method and another method:

restart;
primes:=proc(N) # First way (your method)
local a,b;
b:={};
for a by 2 to N do
if isprime(a)  then b:=`union`(b,{a})
end if;
end do;
return b;
end proc:

primes1:=proc(N) # Second way
local a,b,k;
k:=0:
for a by 2 to N do
if isprime(a)  then k:=k+1; b[k]:=a 
end if;
end do;
return convert(b,set);
end proc:

CodeTools:-Usage(primes(1000000)):

#  memory used=23.16GiB, alloc change=28.49MiB, cpu time=91.33s, real time=96.09s, gc time=3.86s
 

CodeTools:-Usage(primes1(1000000)):

#  memory used=217.08MiB, alloc change=40.00MiB, cpu time=2.28s, real time=2.38s, gc time=171.88ms

Thanks everyone for the helpful comments. I hope this long-standing bug is finally fixed in future Maple releases.

@C_R  If you need to simplify a subexpression (in your example, this is the expression in the numerator under the integral sign), then you can the simplifying procedure apply to the corresponding operand.

restart;
expr := Int(sqrt(-2*alpha^2*x - 2*alpha^2 + 4)/sqrt(-alpha^2 + 1)/(-alpha^2*x - alpha^2 + 2), alpha = 0 .. z);

applyop(e->content(e)*primpart(e), [1,1], expr);

If you open the brackets in the expression  dio1 , you will notice that all of their coefficients for  dio  and  dio1  are the same, except for one.

restart; 
dio := 3*x^2*y^2*z^2-36*x^2*y^2*z-30*x^2*y*z^2-24*x*y^2*z^2+126*x^2*y^2+432*x^2*y*z+93*x^2*z^2+360*x*y^2*z+312*x*y*z^2+66*y^2*z^2-1692*x^2*y-1476*x^2*z-144*x*y^2-5040*x*y*z-1104*x*z^2-1080*y^2*z-948*y*z^2+6120*x^2+21096*x*y+18576*x*z+4554*y^2+16056*y*z+3540*z^2-79848*x-69300*y-61272*z+268452;
dio1:=add(((x-i)*(y-i-1)*(z-i-2))^2, i=1..7,3);
dio2:=expand(dio1);
is(dio-dio1=0);
coeffs(dio);
coeffs(dio2);

This example is fairly easy to solve by hand without Maple. First, we remove the irrationality from the denominator, then simplify the expression under the external root and represent it as a square.

Below are the main step-by-step transformations:

There is no error in Maple 2018.2. Maple simply returns this integral unevaluated:

@janhardo Yes, you're right. The  ~  command isn't needed here.

@sand15  Thank you.

@imparter 
 

restart;
with(plots): with(plottools):
h:=z->1-(delta2/2)*(1 + cos(2*(Pi/L1)*(z - d1 - L1))):
K1:=((4/h(z)^4)-(sin(alpha)/F)-h(z)^2+Nb*h(z)^4):
T1:=Int(K1,z=0..1):
L1:=0.2: F:=10:
d1:=0.2:
T2:=evalf(T1):
T3:=unapply(T2,alpha,Nb,delta2):
T:=textplot([[0.037,4.18,"Color Representation",color=black],[0.04,4.14,"Nb=0.02",color=blue],[0.04,4.1,"Nb=0.03",color=red],[0.04,4.06,"Nb=0.04",color=green]], font=[times,bold,16]):
Lines:=line([0.033,4.14],[0.035,4.14],color=blue,thickness=10),line([0.033,4.1],[0.035,4.1],color=red,thickness=10),line([0.033,4.06],[0.035,4.06],color=green,thickness=10):                     
T11:=textplot([[0.067,4.18,"Line Style"],[0.07,4.14,typeset(alpha=Pi,"/",6)],[0.07,4.1,typeset(alpha=Pi,"/",4),line=longdash],[0.07,4.06,typeset(alpha=Pi,"/",3),line=dash]], font=[times,bold,16]):
Lines1:=line([0.03+0.03,4.14],[0.035+0.03,4.14],linestyle=solid,color=blue),line([0.03+0.03,4.1],[0.035+0.03,4.1],linestyle=longdash,color=red),line([0.03+0.03,4.06],[0.035+0.03,4.06],linestyle=spacedash,color=green):
P:=plot([seq(seq(T3(alpha, Nb, delta2), Nb in [0.1, 0.2,0.3]), alpha in [Pi/6, Pi/4, Pi/3])], delta2 = 0.02.. 0.1,  titlefont = ["ROMAN", 15], labels = ["δ1", "T3"], labeldirections = ["horizontal", "vertical"], labelfont = ["SYMBOL", 11], linestyle = [solid, solid,solid,longdash, longdash,longdash, spacedash,spacedash, spacedash], symbol = [BOX, CROSS, CIRCLE, BOX, CROSS, CIRCLE, BOX, CROSS, CIRCLE], color = [blue,red, green, blue,red, green,blue,red, green], axes = boxed,size = [1000, 1000]):
plots:-display(P,T,T11,Lines,Lines1);

@janhardo  Yes,  minimize  works symbolically, but  the Optimization package - only numerically.