The problem is easily solved by simply enumerating all the options. We get 4 solutions with a minimum distance of 20 :

restart;
pts:=[[0, 0], [1, 1], [2, 5], [4, 2], [5, 3]]:
local D:
A,B,C,D,E:=pts[]:
dist:=(X,Y)->abs(X[1]-Y[1])+abs(X[2]-Y[2]):
P:=combinat:-permute([B,C,D,E]):
k:=0:
for p in P do
k:=k+1; d:=dist(A,p[1])+dist(p[1],p[2])+dist(p[2],p[3])+dist(p[3],p[4])+
dist(p[4],A); L[k]:=[p,d];
od:
L:=sort(convert(L,list),key=(x->x[2])):
L[1..4][];

   L := [[[[1, 1], [2, 5], [5, 3], [4, 2]], 20], [[[1, 1], [4, 2], [5, 3], [2, 5]], 20], [[[2, 5], [5, 3], [4, 2], [1, 1]], 20], [[[4, 2], [5, 3], [2, 5], [1, 1]], 20], [[[1, 1], [2, 5], [4, 2], [5, 3]], 22], [[[1, 1], [5, 3], [4, 2], [2, 5]], 22], [[[2, 5], [4, 2], [5, 3], [1, 1]], 22], [[[5, 3], [4, 2], [2, 5], [1, 1]], 22], [[[1, 1], [4, 2], [2, 5], [5, 3]], 24], [[[1, 1], [5, 3], [2, 5], [4, 2]], 24], [[[4, 2], [2, 5], [5, 3], [1, 1]], 24], [[[5, 3], [2, 5], [4, 2], [1, 1]], 24], [[[2, 5], [1, 1], [4, 2], [5, 3]], 26], [[[2, 5], [1, 1], [5, 3], [4, 2]], 26], [[[4, 2], [5, 3], [1, 1], [2, 5]], 26], [[[5, 3], [4, 2], [1, 1], [2, 5]], 26], [[[2, 5], [5, 3], [1, 1], [4, 2]], 28], [[[4, 2], [1, 1], [2, 5], [5, 3]], 28], [[[4, 2], [1, 1], [5, 3], [2, 5]], 28], [[[5, 3], [2, 5], [1, 1], [4, 2]], 28], [[[2, 5], [4, 2], [1, 1], [5, 3]], 30], [[[4, 2], [2, 5], [1, 1], [5, 3]], 30], [[[5, 3], [1, 1], [2, 5], [4, 2]], 30], [[[5, 3], [1, 1], [4, 2], [2, 5]], 30]]

     [[[1, 1], [2, 5], [5, 3], [4, 2]], 20], [[[1, 1], [4, 2], [5, 3], [2, 5]], 20], [[[2, 5], [5, 3], [4, 2], [1, 1]], 20], [[[4, 2], [5, 3], [2, 5], [1, 1]], 20]

The code below works. 10 frames are allocated for each n :

restart;
WM := n -> plots:-matrixplot(P^n, heights=histogram):
P := Matrix(2$2, [0.8, 0.2, 0.4, 0.6]);
plots:-animate(WM, [n], n=[seq(i$10,i=1..10)]);

restart;
with(LinearAlgebra):
with(GraphTheory):
ts:=time():
g:=Graph(Matrix([[0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0]]));

Eg := add(abs(Eigenvalues(AdjacencyMatrix(g)))):
evalf(%);

DEg := add(abs(Eigenvalues(AllPairsDistance(g)))):
evalf(%);

time()-ts;

 g := Graph 1: an undirected unweighted graph with 31 vertices  and 30 edge(s)
                          37.91134763
                          614.5329657
                             0.922
 

Use the  combinat:-choose  command.

Example:

L:={P1,P2,P3,P4,P5}:
L1:=combinat:-choose(L, 3);
map(mul, L1);

L1 := {{P1, P2, P3}, {P1, P2, P4}, {P1, P2, P5}, {P1, P3, P4}, 

  {P1, P3, P5}, {P1, P4, P5}, {P2, P3, P4}, {P2, P3, P5}, 

  {P2, P4, P5}, {P3, P4, P5}}

 {P1 P2 P3, P1 P2 P4, P1 P2 P5, P1 P3 P4, P1 P3 P5, P1 P4 P5, 

   P2 P3 P4, P2 P3 P5, P2 P4 P5, P3 P4 P5}
 

Use the  tickmarks  option for this.

An example:

restart;
f:=(x,y)->x^2+y^2:
interface(rtablesize=11):
M:=Matrix([seq([seq(f(x,y),x=0..5,0.5)],y=0..5,0.5)]);
plots:-matrixplot(M, heights = histogram, colorscheme = ["Blue", "Green", "Yellow", "Red"], axes=normal, tickmarks=[[seq(i=0.5*i-0.75,i=1.5..11.5,1)],[seq(i=0.5*i-0.75,i=1.5..11.5,1)],default], labels=[x,y,"f(x,y)"], orientation=[-75,75]);

Looks like a bug. Below is a workaround:

plots:-display(
	plots:-implicitplot(`if`(y<=1-x^2,-x^2+y,undefined), x = 0 .. 2, y = -2 .. 1, color=red),
	plot(1-x^2, x=0..2, color=blue), scaling=constrained, view=[0..2,-2..1]
);

In fact, your expression is a first degree polynomial with respect to  exp(x) . Its coefficients should be simplified but not factored:

restart;
expr:= (cos(x)^2+sin(x)^2)+5+(1+x+x^2+x^3)*(cos(x)^2+sin(x)^2)*exp(x);
simplify~([coeffs(expr,exp(x),'t')]);
add(%*~t);

Another reliable way to find all complex roots (in particular, real roots) of some analytic function in the specified intervals is to use RootFinding:-Analytic command:

restart:
line := x/100 - 1/2:
wave := cos(x/5) * sin(x/2): ## -1 <= wave <= 1
eq:= line - wave:
L:=sort([RootFinding:-Analytic(eq, re=-50..150, im=-1..1)]);
nops(L);
plot(eq, x=-50..150, -0.05..0.05, color=red, size=[1000,400]);

The line

((max-min)/abs(Mean))(S);

is incorrect. Should be

(max-min)(S)/abs(Mean(S)); 

The result is the same as 

((max-min)/(abs@Mean))(S);

You have one equation with three unknowns. Such equations usually have infinitely many solutions. The plot clearly shows that your equation is satisfied by the points lying on  some surfaces:

plots:-implicitplot3d((x,y,z)->`if`(y>x and z>y,cos(2*Pi*(x+y-2*z))+cos(2*Pi*(y+z-2*x))+cos(2*Pi*(z+x-2*y)),undefined),  0..1, 0..1, 0..1, style=surface, grid=[50,50,50], axes=normal);

restart;
h:=k-(k-1)*x: lambda:=k/(1+k):
DE:=diff(p(x),x)=6/h^2-12/h^3*lambda;
BC:=p(0)=0, p(1)=0;
simplify(dsolve({DE, BC}, p(x)));
is(eval(p(x),%)=6*(h-1)*(h-k)/h^2/(1-k^2));

Or

restart;
h:=k-(k-1)*x: lambda:=k/(1+k):
DE:=diff(p(x),x)=6/h^2-12/h^3*lambda;
BC:=p(0)=0, p(1)=0;
simplify(dsolve({DE, BC}, p(x)));
p:=unapply(eval(p(x),%), x);
h:='h':
eval(p(x),x=solve(h=k-(k-1)*x,x));
simplify(%);

int(u*v, x=-1..0, y=-1..1) + int(u*v, x=0..1, y=-1..1);

expr:=sin(x)+cos(x)*exp(x)+x*tan(x)*f(x)+x+y;
indets(expr, function) union {x};

                         {x, cos(x), exp(x), f(x), sin(x), tan(x)}

restart;
plot3d([[r*cos(phi),r*sin(phi),r],[r*cos(phi),r*sin(phi),2-sqrt(4-r^2)]], r=0..2, phi=0..2*Pi, scaling=constrained);

Another option:

plot3d([[r*cos(phi),r*sin(phi),r],[r*cos(phi),r*sin(phi),2+sqrt(4-r^2)]], r=0..2, phi=0..2*Pi, scaling=constrained);

@flesx  R is your initial region.

restart;
R:=plot([seq(k,k=0..1,0.1),seq([k,y,y=0..1],k=0..1,0.1)],x=0..1, color=blue, scaling=constrained):
z:=x+I*y:
F:=plottools:-transform(unapply(evalc([(Re,Im)(exp(2*z))]), x,y)):
plots:-display(<R | F(R)>);