@Carl Love  Unfortunately I have the right to put only 1 vote up.

@Carl Love  for a very useful remark. I fixed this defect in code.

Please provide explicitly your system  V . Probably needed extra options in  solve  command.

@Carl Love  Of course, in the original example, there is no problem, because determinant of the system at any  b  is not equal to 0, and it has a unique solution for any b. In other cases (if first we use solve) there may be a wrong solution.

Example:

restart;

solve({b*x+y=1, x+b*y=1}, {x,y});

eval(%, b=1);

The solution is erroneous because for b = 1, there are infinitely many solutions.

The correct solution

solve(eval({b*x+y=1,x+b*y=1}, b=1), {x,y});

                   {x = 1-y, y = y}

Another example:

eval(solve(sqrt(x+b)>x, x), b=1);

Warning, solutions may have been lost

The correct solution:

solve(eval(sqrt(x+b)>x, b=1), x);

               RealRange(-1,Open(1/2*5^(1/2)+1/2))

Thus, variant  solve(eval...)  is more reliable than  eval(solve...)

@Carl Love  for clarifying the situation.

@Carl Love 

Q:=evalf[20](P);

L:=[fsolve(Q, complex)];

seq(abs(subs(x=t, P)), t=L);

What does  HFloat(undefined)  mean?

@acer  Thank you for your helpful advice. Unfortunately, the situation is complicated by the fact that in the case where the points lie on parallel faces, pagan's way and accordingly, my procedure may give incorrect results.

Here is an example where the shortest path passes through four faces:

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

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]);