I don't know any other way to solve the problem than to do it using the tools of the plots and plottools packages:

restart; 
with(plots): with(plottools): 
P := plot(sin(x), x = -3 .. 3, colour = [blue], style = pointline, symbol = solidcircle, numpoints = 20): 
p1 := line([1.5, -1], [1.9, -1], colour = [blue]): 
p2 := pointplot([1.7, -1], symbol = solidcircle, color = blue, symbolsize = 10): 
p3 := curve([[1.4, -.9], [2.65, -.9], [2.65, -1.1], [1.4, -1.1], [1.4, -.9]]): 
p4 := textplot([2.3, -1, sin(x)]): 
plots:-display(P, p1, p2, p3, p4);

Download Plot1_new.mw

The problem can be easily reduced to solving a system of 2 equations with 2 unknowns. But for some reason, there is no correct answer among the suggested ones (possibly a typo):

restart;
d:=(X,Y)->sqrt((X[1]-Y[1])^2+(X[2]-Y[2])^2):
A:=[7,6]: B:=[3,4]: P:=[x,0]: Q:=[0,y]:
solve({d(P,A)=d(P,B),d(Q,A)=d(Q,B)});
assign(%);
d(P,Q); # The answer

@vv You only check six-digit numbers with distinct digits. But how do we know that numbers with fewer digits and/or repeating digits cannot have this property?

The following code checks for any numbers that have no more than six digits:

restart;
L:=[$2..6]:
for k from 2 to 6 do
for n from 10^(k-1) to floor((10^k-1)/6) do
L1:=convert(n,base,10);
S:=combinat:-permute(L1);
P:=n*L;
L2:=convert~(P,base,10);
if andmap(t->t in S, L2) then print(n) fi;
od: od:

                                                    142857

If I understood your task correctly, here is my attempt:

restart;
P:=proc(x1,p2,t)
local v, d12, v12, r, x1_t, C, c, d1, s, T;
uses plots,plottools;
v:=convert(p2-x1,Vector);
d12:=sqrt(v.v);
v12:=v/d12;
x1_t:=convert(x1,Vector)+t*v12;
r:=d12-t;
d1:=pointplot(convert(x1_t,list),color=red,symbol=solidcircle,symbolsize=12);
C:=circle(x1,d12, linestyle=3, color=blue);
c:=circle(convert(x1_t,list),r,color=blue);
s:=line(x1,p2, thickness=2);
T:=textplot([[x1[],"x1",align=[left,above]], [p2[],"p2", align=[right,above]]], font=[times,15]):
display(d1,c,s,T,C, scaling=constrained);
end proc:

plots:-animate(P,[[1,1],[2,3],t], t=0..1, frames=60);

The problem is solved by direct calculation using the shoelace formula. See  https://en.wikipedia.org/wiki/Shoelace_formula

restart;
local D:
A:=[0,0]: B:=[b,c]: C:=[a+b,c]: D:=[a,0]:
f:=(A,B)->sort(expand((y-A[2])*(B[1]-A[1])-(x-A[1])*(B[2]-A[2]))):
g:=(a,b)->eval([x,y],solve({a, b},{x,y})):
mAB:=(A+B)/2: mBC:=(B+C)/2: mCD:=(C+D)/2: mDA:=(D+A)/2:
A1,A2,B1,B2,C1,C2,D1,D2:=map(t->f(t[]),[[A,mCD],[A,mBC],[B,mDA],[B,mCD],[C,mAB],[C,mDA],[D,mBC],[D,mAB]])[]:
L:=[[A2,B1],[A2,C1],[B2,C1],[B2,D1],[C2,D1],[C2,A1],[D2,A1],[D2,B1]]:
assign(seq(p[i]=g(L[i][]), i=1..8));
p[9]:=p[1]:
S:=expand(1/2*add(p[i][1]*p[i+1][2]-p[i][2]*p[i+1][1], i=1..8)):
is(abs(S)=abs(a*c/6));

 #  true

Explanations of the code:
A, B, C, D are the vertices of an arbitrary parallelogram (a<>0, c<>0); mAB, mBC, mCD, mDA are the midpoints of the corresponding sides, f is the procedure for finding a straight line through 2 points, g is the procedure for finding the intersection point of two straight lines, p[i] (i=1..8) are the vertices of the octagon, S is the area of ​​the octagon. The calculations use only formulas from the Maple core, so they are suitable for any version of Maple (even very old ones).

Visualization - the names of the objects match their names in the code above:

                           

Edited.                       

Octagon.mw

Since all the numbers in the problem are multiples of 5, the idea immediately came to me to work with a regular pentagon (roses - red solid circles and rows - blue segments):

restart;
with(plots): with(plottools):
Line:=(A,B)->sort(expand((y-A[2])*(B[1]-A[1])-(x-A[1])*(B[2]-A[2]))):
assign(seq(A[i]=[cos(Pi/2+2*Pi*(i-1)/5),sin(Pi/2+2*Pi*(i-1)/5)],i=1..5));
B[1]:=eval([x,y],solve({Line(A[3],A[5]),x=0})):
assign(seq(B[i]=[abs(B[1][2])*cos(-Pi/2+2*Pi*(i-1)/5),abs(B[1][2])*sin(-Pi/2+2*Pi*(i-1)/5)],i=2..5)):
E:=(B[3]+B[4])/2:
C:=eval([x,y],solve({Line(B[3],B[5]),Line(B[2],B[4])})):
M:=eval([x,y],solve({Line(B[1],B[4]),Line(B[2],B[5])})):
N:=eval([x,y],solve({Line(B[1],B[3]),Line(B[2],B[5])})):
K:=(M+N)/2:
P1:=pointplot([A[1],E,C,K,B[1]], color=red, symbol=solidcircle, symbolsize=15):
P2:=display(line(A[1],A[3]),line(A[1],B[1]),line(A[1],A[4]),line(B[2],B[5]), color=blue):
display(seq(rotate(P1,2*Pi*k/5),k=1..5),seq(rotate(P2,2*Pi*k/5),k=1..5), axes=none, size=[650,650], scaling=constrained);

Edited.

restart;
M:=<sqrt(5),0>: C:=<2*sqrt(5),0>: MB:=y=2*(x-sqrt(5)):
Circle:=x^2+y^2=(2*sqrt(5))^2:
solve({MB,Circle}):
E:=eval(<x,y>,%[2]):
EM:=M-E: EC:=C-E:
Angle_CEM=arccos(EM.EC/sqrt(EM.EM)/sqrt(EC.EC));

In cases where the  solve command fails, it is useful to first plot the graphs to visually verify the roots. Then we find the roots numerically using the  fsolve  command. Finally, we try to find the roots symbolically using the  identify  command:

Download roots.mw

a:=Matrix(2,2,{(1, 1) = <1.0,2.0>, (1, 2) = 3.0, (2, 1) = 4.0, (2, 2) =<5.0,6.0>});
a[1,1][1];
a[2,2][2];

We will assume that the parameters are not equal to  0  simultaneously, otherwise we have a trivial case  x=0 . Factoring, we see that for any values ​​of the parameters, the equation has the root  n=1 . The second factor is a third-degree polynomial with respect to n  with the leading coefficient different from 0, so there will be another real root:

Thus, for any values ​​of the parameters, the equation will have at least  2  real roots.

Download quartic_equation_in_n_new.mw

Let's consider the vertices of a regular hexagon with a side of 1, the center of which lies at the origin. For definiteness, we will assume that the point in the center is red. Since there are still 6 points left, there are a total of 2^6=64 coloring options (see the figure below). Of these 64 options, only 3 options do not have an equilateral triangle with equally colored vertices. Let's take one of these options (the next figure) and add 3 more points to the right, forming 2 equilateral triangles with a side of 1 with the right red point. If the rightmost point is black, then we immediately get an equilateral triangle with black vertices with a side  sqrt(3). If this point is red, then with any choice of colors for the remaining 2 points we also get an equilateral triangle with vertices of the same color.

restart;
with(plots):
L:=combinat:-permute([red$6, black$6], 6):
P:=seq(pointplot([[0,0],seq([cos(2*Pi*k/6),sin(2*Pi*k/6)], k=0..5)], color=[red,p[]],symbol=solidcircle, symbolsize=40, scaling=constrained, axes=none, size=[100,100]), p=L):
display(Matrix(8,8,[P]));
pointplot([[0,0],seq([cos(2*Pi*k/6),sin(2*Pi*k/6)], k=0..5),[3/2,sqrt(3)/2],[2,0],[3/2,-sqrt(3)/2]], color=[red,L[28][],"LightGrey","LightGrey","LightGrey"], scaling=constrained,symbol=solidcircle, symbolsize=20, size=[600,600]);

Edit. You can get by with just one additional point above the top side of the hexagon, as  @Alfred_F  already wrote about. Regardless of the choice of its color, we get an equilateral triangle with vertices of the same color:

Look at the following 2 examples:

cat(alpha,beta); # not OK
`αβ`; # OK
cat(d,alpha); # not OK
`dα`; # OK

For your example

`#3  Q24` (`α`);

restart;
Proc := proc(m)
local xA, yA, xB, yB, xC, yC, xJ, yJ, tx, Oo, c1, r, eqBJ, eq1, sol, O;
uses plots, geometry;
_EnvHorizontalName := 'x'; _EnvVerticalName := 'y';
xJ := 5; yJ := 1; point(A, 2, 4); point(J, xJ, yJ); point(Oo, 0, 0);
r := 3; circle(c1,[Oo, r]);
eqBJ := y = m*(x - xJ) + yJ; line(BJ, eqBJ, [x, y]);
eq1 := x^2 + y^2 = r^2; sol := solve({eqBJ, eq1}, {x, y}, explicit);
xB := subs(sol[1], x); yB := subs(sol[1], y);
point(B, xB, yB); xC := subs(sol[2], x); yC := subs(sol[2], y);
point(C, xC, yC); circle(c2, [A, B, C],centername=O); line(AB, [A, B]); line(AC, [A, C]); line(BC, [B, C]);
eqBJ := y = m*(x - xJ) + yJ; line(BJ, eqBJ, [x, y]);
eq1 := x^2 + y^2 = r^2; sol := solve({eqBJ, eq1}, {x, y},explicit);
xB := subs(sol[1], x); yB := subs(sol[1], y); point(B, xB, yB);
xC := subs(sol[2], x); yC := subs(sol[2], y); point(C, xC, yC);
line(AB, [A, B]); line(AC, [A, C]);
draw([BC(color = black), c2(color = magenta), A(color = blue, symbol = solidcircle, symbolsize = 16),
B(color = red, symbol = solidcircle, symbolsize = 16), C(color = red, symbol = solidcircle, symbolsize = 16),
J(color = red, symbol = solidcircle, symbolsize = 16),c1(color=blue)],axes = normal, printtext = true);
end proc:

plots:-animate(Proc,[m],m=-0.3..1, frames=50);  # Animation

plots:-display(Proc(-0.3),Proc(0),Proc(1)); # 3 circles

It is clearly seen that these circles have only one common point  A .

Proc.mw

The formula  cos(theta) = sqrt(1 - sin(theta)^2)   is incorrect because the cosine can be negative, but the root on the right is always non-negative. The right  substitution  is  cos(theta) = sin(theta + Pi/2)  or  cos(theta) = sin(Pi/2  -  theta ) .

restart;
eqs:=[a=b,c=d,e=f]:
map(`+`, eqs, 2);
map(`-`, eqs, 2);
map(`/`, eqs, 2);
map(t->t*~2, eqs);
map(t->t^2, eqs);