I don't know why Student:-Calculus1:-VolumeOfRevolution command distorts the colors. You may get better results if you use the  plot3d  command to draw surfaces. I also selected the edges of the cylinder (2 circles) using the  plots:-spacecurve  command:

restart;
A:=plot3d([x,4*cos(t),4*sin(t)], t=0..2*Pi, x=1..2, style=surface, color="Blue"):
B:=plot3d([[1,r*cos(t),r*sin(t)],[2,r*cos(t),r*sin(t)]], t=0..2*Pi, r=0..4, style=surface, color="Blue"):
C:=plots:-spacecurve([[1,4*cos(t),4*sin(t)],[2,4*cos(t),4*sin(t)]], t=0..2*Pi, color=black, thickness=2):
plots:-display(A,B,C, transparency=0.6, scaling=constrained, view = [0 .. 3, -5 .. 5, -5 .. 5], axes=normal, orientation = [-94,-18,11], labels = [x, y, z]); 

See the following toy example:

A:=<1,2; ``,3>;
map(t->`if`(t=``,``,f(t)), A);

Look at these 2 graphs. We see that on a part of the range  0..2*Pi , the function takes negative values. Therefore, the square root on these intervals will take on complex values.

plot(c^2-a^2*cos(theta)^2, theta = 0 .. 2*Pi); 
plot(b^2-(a*sin(theta)+sqrt(c^2-a^2*cos(theta)^2))^2, theta = 0 .. 2*Pi);

You must use the  solve  command for this:

restart;
eq:=5*a^2 + 2*b - 7 = 0;
a:=1;
b=solve(eq, b);

This can be done more simply and more automatically as in the example:

restart;
with(plots):
with(plottools):
p1:= cuboid( [0,0,0], [1,1,1], color="LightBlue", transparency=0.5): # The cube
p2:=plot3d(3/2-x-y, x=-1..2, y=-1..2, style=surface, color=red): # The plane
display([p1, p2], view=[(0..1)$3]);

I considered a more interesting example in which the intersection will be a regular hexagon. For your example in the code replace 3/2 - x - y  with  1 - x .

The problem is surprisingly easy to solve. Leave the first column unchanged. For the remaining columns, sequentially in the cycle, we find the desired permutation:

restart;
L1:=[blue$3,brown$2,yellow$2,red,black]:
L2:=[blue$3,purple$2,pink$2,red,black]:
L3:=[green$3,brown$2,pink$2,red,black]:
L4:=[orange$3,purple$2,yellow$2,red,black]:
L5:=[green$3,brown$2,yellow$2,red,black]:
L6:=[orange$3,purple$2,pink$2,red,black]:
L7:=[blue$3,brown$2,yellow$2,red,black]:
assign(seq(P||i=combinat:-permute(L||i), i=2..7)):

OneStep:=proc(L,P)
local p;
for p in P do
if `and`(seq(not (p[i] in convert(L[i],set)), i=1..9)) then
return [seq([op(L[i]),p[i]],i=1..9)] fi;
od;
end proc:

L:=L1:
for n from 1 to 6 do
L:=OneStep(L,P||(n+1));
od:

for n from 1 to 9 do
R[n]:=plots:-pointplot([seq([i,10-n],i=1..7)], color=convert(L[n],list), symbol=solidcircle, symbolsize=40);
od:
plots:-display(seq(R[n], n=1..9), axes=none, size=[500,700], scaling=constrained); 

I don't know of a purely automatic way to solve the problem. The procedure below returns a set of self-intersection points. But it uses an optional parameter  n  (default is 50). The number  n  means the number of intervals for applying the  fsolve  command. We select the number  n  in such a way that only 1 self-intersection point can fall into one interval. I think that for most examples  n=50  will be enough (in the example below,  n=30  was enough). But probably there are complex examples where a larger number of  n  is required. Therefore, it is recommended to make a plot before applying the procedure in order to estimate the value of  n .

restart;
Self_Intersection_Points:=proc(L::list,R::range,n::posint:=50)
local x, y, a, b, h, k, m, Sol, S;
x:=unapply(L[1],t); y:=unapply(L[2],t);
a:=lhs(R); b:=rhs(R);
h:=(b-a)/n;
m:=0;
for k from 2 to n do
Sol:=fsolve({x(t)=x(s),y(t)=y(s)}, {s=a+h*(k-1)..a+h*k,t=a..a+h*(k-1)});
if type(Sol,set) then m:=m+1; S[m]:=eval(s,Sol) fi;
od;
convert(S,set);
end proc:

Example of use:

x:=t->4*cos(t/2)+2*cos(2*t)+cos(4*t): y:=t->4*sin(t/2)+2*sin(2*t)+sin(4*t):
L:=Self_Intersection_Points([x(t),y(t)],0..4*Pi, 30);
A:=plot([x(t),y(t), t=0..4*Pi], color=blue):
B:=plot(map(t->[x(t),y(t)],L), style=point, symbol=solidcircle, color=red, symbolsize=10):
plots:-display(A,B, size=[500,500]);

Your decision is correct. You can test each of the solutions using the  odetest  command:

restart;
Dgl_1 := diff(y(x), x) = (1 + y(x)^2)/(x*y(x)); 
P := dsolve(Dgl_1, y(x));
odetest(P[1], Dgl_1);
odetest(P[2], Dgl_1); 

restart;
f(x)=piecewise(x<>0,5,x=0,0);
# or
f(x)=piecewise(x<>0,5,0);

See the post  https://www.mapleprimes.com/posts/202222-Contour-Curves-With-Labels

restart;
P:=3*x^2 + a*x*y + b*y^2:
var:=[x,y]:
S1:=[coeffs(P,var,t)];
S2:=[t];
[seq([S1[i],[degree(S2[i],x),degree(S2[i],y)]], i=1..nops(P))];
# Or (shorter)
zip((u,v)->[u,[degree(v,x),degree(v,y)]], S1,S2);

                                       S1 := [3, a, b]
                                  S2 := [x^2, x*y, y^2]
                         [[3, [2, 0]], [a, [1, 1]], [b, [0, 2]]]
                         [[3, [2, 0]], [a, [1, 1]], [b, [0, 2]]]

The  Student:-Calculus1:-ShowSolution  command does not work with complex expressions. But you can (using the  evalc  command) first split the integral into real and complex parts (assuming  p  is a real number), and then use the  Student:-Calculus1:-ShowSolution  command separately for each part:

Download ShowSolution.mw

If you use the plot command to plot points, then add the  style=point  option. Otherwise, Maple connects consecutive points with line segments (if several consecutive points lie on the same straight line, then they will not be visible on the plot):

restart;
Points:=[[10, 2], [11, 2], [12, 2], [13, 2], [14, 2], [15, 2], [16, 2], [17, 3], [18, 2], [19, 2], [20, 2], [21, 3], [22, 3], [23, 3], [24, 3], [25, 2], [26, 3], [27, 3], [28, 3], [29, 3], [30, 3], [31, 4], [32, 3], [33, 3], [34, 3], [35, 3], [36, 3], [37, 4], [38, 4], [39, 4], [40, 4], [41, 3], [42, 3], [43, 4], [44, 4], [45, 4], [46, 4], [47, 4], [48, 4], [49, 4], [50, 4], [51, 4], [52, 4], [53, 4], [54, 4], [55, 4], [56, 4], [57, 5], [58, 5], [59, 5], [60, 5], [61, 4], [62, 4], [63, 4], [64, 4], [65, 5], [66, 5], [67, 5], [68, 5], [69, 5], [70, 5], [71, 5], [72, 4], [73, 5], [74, 5], [75, 5], [76, 5], [77, 5], [78, 5], [79, 5], [80, 5], [81, 5], [82, 6], [83, 6], [84, 6], [85, 5], [86, 5], [87, 5], [88, 5], [89, 5], [90, 5], [91, 6], [92, 6], [93, 6], [94, 6], [95, 6], [96, 6], [97, 6], [98, 5], [99, 5], [100, 5]]:
nops(Points);
plot(Points, style=point, color=red, symbol=solidcircle, symbolsize=5, size=[1000,400]);

All 91 points were built successfully.

It is more natural to use 2 indices, where the first index is the degree of the first variable, and the second index is the degree of the second variable:

P:=(x,y,N)->add(add(c[i,j]*x^i*y^j, i=0..N-j),j=0..N):

sort(P(x,y,5),[x,y]);  # Example of use

See help  ?ScientificConstants  and  ?initialconstants .