I found an error in my solution (Pagan's approach also wrong). If two points lie on the adjacent faces, the shortest path is not necessarily passes only through these faces (it may pass through the 3 faces).

Here is the example:

A:=plottools[cuboid]([0,0,0], [1,1,1],style=wireframe, color=black):
B:=plottools[curve]([[1,1/2,1/10],[1,3/5,0],[3/5,1,0],[1/2,1,1/10]], color=red, thickness=3):
C:=plots[textplot3d]([[1,1/2,1/10,"A"],[1/2,1,1/10,"B"]], color=black, align=left, font=[TIMES,ROMAN,18]):
plots[display](A,B,C);

I will think how to fix this situation.

@Preben Alsholm  I just wanted to remind that the conversion a list to a set - a convenient way of getting rid of duplicate elements. Simultaneously, Maple automatically sorts the remaining elements (lexicographically).  If the original order is irrelevant, you can save this set without turning it back to the list:

u:=[x,y,z]; v:=[a,b,x];

{op(u), op(v)};

              {a,b,x,y,z}

@Markiyan Hirnyk  And how to use it to generate a specific random matrix?

@Markiyan Hirnyk  I did not understand that your command returns.

M := Matrix(3, 5,  (i, j) -> Statistics:-RandomVariable(Normal(0, 1)));

@acer  Yes, your method is much faster for large matrices (vote up). I wonder what the reason for the slower work  of  RandomTools[Generate]

@Kanellopoulos  The procedure with 2 arguments:

L[1]:=[x1,x2,x3,x4]:
L[2]:=[y1,y2,y3,y4]:

K:=(j::posint, i::integer) -> `if`(irem(i,4)<>0, L[j][irem(i,4)], L[j][4]):

Examples:

K(1,4), K(2,5), K(2,0);

                 x4, y1, y4

Addition: If you want to make this working also for negative  i , replace  irem(i,4)  by  i mod 4

@erik10

PlotFootball := proc(L::list)

local FootballFaces;

uses geom3d, plottools, plots;

  TruncatedIcosahedron(football, point(C, 0, 0, 0), 1);

  FootballFaces := seq(polygon(faces(football)[i], color = `if`(nops(faces(football)[i])=5,L[1],L[2]), linestyle=solid, tickness=3), i = 1 .. 32);

display(FootballFaces, axes = none);

end proc: 

Example:

PlotFootball([red,yellow]);

Maple plots  graphs of expressions (rather than functions) in their domains. The domain of expression - it is the set of all the values of a variable in which the value of the expression is a real number. When operations on the expressions their domains can be changed.

At a strict definition of the function  its domain should be included in the definition. Maple is not obliged to do this for you (it just does not know this). So the workaround for your plotting is

plot([sqrt(x+2*sqrt(x-1))+sqrt(x-2*sqrt(x-1)), sqrt(x-2*sqrt(x-1)), sqrt(x+2*sqrt(x-1))], x = 1 .. 5, color=[red,blue,green], thickness=2,numpoints=5000, scaling=constrained, view=[-0.5..5.5, -0.5..4.5]);

@maple2015  Because Maple simplifies  (sqrt(x-1)+1) + (1-sqrt(x-1)) = 2  and forgets about the domain x>=1 . Of course it was a mistake Maple.

Addition: The similar example

is(x/x=1);

    true

Maple automatically simplifies  x/x  into  1  and  forgets that  x  should be not 0

@Markiyan Hirnyk   Your recent comments like trolling, so I stop this discussion.

@Markiyan Hirnyk   Because  x= 0.02612014634  does not belong to  0.01338216460 < x and x <= 0.02612014630

@Markiyan Hirnyk  For  i=14

fsolve(p[14][2] = VP[14]);

                0.02612014634

Compare with the range  0.01338216460 < x and x <= 0.02612014630

For  i=2

fsolve(p[2][2] = VP[2]);

       0.4793595141

The corresponding range  0.3166653569 < x and x <= 0.4793595141

@acer  Thank you very much for the detailed analysis of the situation with floats.

@Markiyan Hirnyk  Your example does not contradict anything. By default  Digits=10 so

is(4.=4.+10^(-20));

              true

Or use quotes:

plot( 'myproc(300,V)', V=0.1..0.9 );

Maximize(
'myproc(300,v)',
v=0.00001..0.99999
);