@Rouben Rostamian   Here is a shorter version of the same plotting without  plots  package:

x := [0., .20, .40, .60, .80, 1.0]:
y1 := [0., .20, .40, .60, .80, 1.0]:
y2 := [0., .64, .96, .96, .64, 0.]:

plot([`[]`~(x,y1),`[]`~(x,y2)]);

I think I misunderstood the meaning of the original question in my answer.

See the same question

@toandhsp 

restart:

k:=combine((sin(x+Pi/3))^2);

sol:=solve(k=1/4, x, AllSolutions = true, explicit);

sol1, sol2:=subs({_B1=0,_Z1=n},sol), subs({_B1=1,_Z1=n},sol);

@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]