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 .
 

seq(2*(n-1)/(2*n-1), n=1..10);
seq((7*n-5)/(7*n-4), n=1..10);
seq((2*n-1)^2/(2*n)^2, n=1..10);
seq((n^2-1)/n^2, n=1..10);

Download SystemRecursive_new.mw

Download System_Recursive_Equations_new.mw

restart;	
solve({2*(k+1)*X(k+1)+(k+1)*Y(k+1)-X(k)-Y(k)+(1)/k!,(k+1)*X(k+1)+(k+1)*Y(k+1)+2*X(k)+Y(k)+(1)/k!},{X(k+1),Y(k+1)});
sol:=subs(k=k-1,eval([X(k+1),Y(k+1)],%));

Sol:=proc(n)
option remember;
if n=0 then return [2,1] else eval(sol,[X(k-1)=Sol(n-1)[1],Y(k-1)=Sol(n-1)[2],k=n]) fi;
end proc:

seq(Sol(k), k=0..5); # Example

The output in the example above  was corrected.

The  evalc  command performs operations on complex expressions, assuming all parameters are real numbers. The  normal  command brings the sum (difference) of fractions to a common denominator:

restart;
Zin := Rin + omega*L*I + (omega*L)^2/(RL + omega*L*I);
evalc(Zin);
normal(Re(%))+normal(Im(%))*I assuming real;