The procedure  Dist  checks the points  P1  and  P2  whether they belong to the surface of the cube  x=0..1, y=0..1, z=0..1, and then  (following pagan's geometric method), finds the distance between them.

Dist:=proc(P1::list(numeric),P2::list(numeric))

local IsOnCube, d1, d2, d3, d4;

uses geom3d;

IsOnCube:=proc(x,y,z)

if (x>=0 and x<=1) and (y>=0 and y<=1) and (z>=0 and z<=1) and ((x=0 or x=1) or (y=0 or y=1) or (z=0 or z=1)) then return true else false fi;

end proc;

if IsOnCube(op(P1))=false or IsOnCube(op(P2))=false then error "P1 or P2 don't lie on the cube" fi;

point(Q1,P1);  point(Q2,P2);

if ycoord(Q1)<>0 then

  if xcoord(Q1)=0 then rotation(Q1,Q1,Pi/2,line(l,[1/2,1/2,t],t));     rotation(Q2,Q2,Pi/2,line(l,[1/2,1/2,t],t)) else

  if ycoord(Q1)=1 then rotation(Q1,Q1,Pi,line(l,[1/2,1/2,t],t));     rotation(Q2,Q2,Pi,line(l,[1/2,1/2,t],t)) else

  if xcoord(Q1)=1 then rotation(Q1,Q1,-Pi/2,line(l,[1/2,1/2,t],t));     rotation(Q2,Q2,-Pi/2,line(l,[1/2,1/2,t],t)) else

  if zcoord(Q1)=0 then rotation(Q1,Q1,-Pi/2,line(l,[t,1/2,1/2],t));     rotation(Q2,Q2,-Pi/2,line(l,[t,1/2,1/2],t)) else

  if zcoord(Q1)=1 then rotation(Q1,Q1,Pi/2,line(l,[t,1/2,1/2],t));     rotation(Q2,Q2,Pi/2,line(l,[t,1/2,1/2],t)) fi;fi;fi;fi;fi;

fi;

if ycoord(Q2)=0 then return distance(Q1,Q2) else

  if xcoord(Q2)=0 then return distance(Q1,rotation(Q2,Q2,Pi/2,line(l,[0,0,t],t))) else

  if zcoord(Q2)=0 then return distance(Q1,rotation(Q2,Q2,-Pi/2,line(l,[t,0,0],t))) else

  if xcoord(Q2)=1 then return distance(Q1,rotation(Q2,Q2,-Pi/2,line(l,[1,0,t],t))) else

  if zcoord(Q2)=1 then return distance(Q1,rotation(Q2,Q2,Pi/2,line(l,[t,0,1],t))) fi;fi;fi;fi;fi;

if ycoord(Q2)=1 then

    d1:=distance(Q1,rotation(Q2,rotation(Q2,Q2,Pi/2,line(l,[0,1,t],t)),Pi/2,line(l,[0,0,t],t)));

    d2:=distance(Q1,rotation(Q2,rotation(Q2,Q2,-Pi/2,line(l,[t,1,0],t)),-Pi/2,line(l,[t,0,0],t)));

    d3:=distance(Q1,rotation(Q2,rotation(Q2,Q2,-Pi/2,line(l,[1,1,t],t)),-Pi/2,line(l,[0,0,t],t)));

    d4:=distance(Q1,rotation(Q2,rotation(Q2,Q2,Pi/2,line(l,[t,1,1],t)),Pi/2,line(l,[t,0,0],t))); 

return min(d1,d2,d3,d4);

fi;

end proc:

Examples:

Dist([1/2, 1/2, 0], [1, 1, 1]);

                    (1/2)*sqrt(10)

Dist([1/2, 1/2, 0], [1/2, 1/2, 1]);

                                2

Dist([0, 1, 0], [1, 1, 1]);

                           sqrt(2)

Distance_between_points.mw

If you want to have the resulting list has been sorted, then

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

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

Matrix(5, RandomTools[Generate](distribution(Normal(10,3)), makeproc=true));

L:=[x1,x2,x3,x4]:

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

Examples:

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

                  x4, x1, x2, x3, x4, x1

I have simplified and slightly changed your procedure: tr  removed and added the list  L  of colors for faces:

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 = L[i]), i = 1 .. 32);

display(FootballFaces, axes = none);

end proc:

Example with 4 colors:

PlotFootball([yellow$8, red$8, blue$8, green$8]);

To maintain order, use square brackets instead of curly ones:

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

Green is  sqrt(x-2*sqrt(x-1)) = sqrt(x-1)+1

Blue is  sqrt(x-2*sqrt(x-1)) = abs(sqrt(x-1)-1) = piecewise(x<=2, 1-sqrt(x-1), x>2, sqrt(x-1)-1)

Red = blue + green

I wrote on that Russian forum that the essence of the problem is that the function  g(x)  is the broken line, defined by the points  seq([VE[i], VP[i]], i=1..15)  .  Therefore, the solutions of the equation  g(x)=VP[i]  for each interval  VE[i+1] < x <= VE[i]  coincides with the right boundary of the interval ,  i.e. with  VE[i] . Due to roundoff errors (as Carl wrote) these numbers can be in the following range, so in the range  VE[i+1] < x <= VE[i]  there no solutions.

It is natural to generalize the problem to find all the positions of an element  in a Matrix. The simple procedure  Positions  makes  it. The procedure returns the sequence of number of occurrences of the element  a  in the matrix  A  and the list of all positions of this element.

Positions:=proc(a, A::Matrix)

local k, m, n, i, j, L;

k:=0;

m, n:=op(1,A);

for i to m do

for j to n do

if A[i,j]=a then k:=k+1; L[k]:=[i,j] fi;

od; od;

if k=0 then k,[] else k,convert(L,list) fi;

end proc:

Example of use:

A:=LinearAlgebra[RandomMatrix](5,6, generator=1..9);

 Positions(4, A);

@Rouben Rostamian  Excellent solution  (vote up). You wrote about presentation as a Maplet. I see no difference in principle between this and a procedure. In the latter the user can also specify any valid parameters.

Formal parameters in the procedure Pursuit :  d  is the initial distance between the pirate ship and the merchant ship,  v1  and  v2  are their velocities, theta0  is the polar angle the motion of the merchant ship.

The code of the procedure:

Pursuit := proc (d, v1, v2, theta0)

local b, r, x, y, t0, theta, x1, y1, x2, y2, T, A, B, P1, P2, A1, B1;

uses plots;

b := solve(sqrt(b^2+1)/b = v1/v2);

r := unapply(d*v2*exp(b*theta)/(v1+v2), theta);

x := unapply(r(theta)*cos(theta), theta);

y := unapply(r(theta)*sin(theta), theta);

t0 := d/(v1+v2);

dsolve({(diff(x(theta(t)), t))^2+(diff(y(theta(t)), t))^2 = v1^2, theta(t0) = 0});

theta := unapply(op(select(x->is(0 < eval(diff(x, t), t = t0)), [rhs(%[1]), rhs(%[2])])), t);

x1 := unapply(piecewise(0 <= t and t < t0, d-v1*t, x(theta(t))), t);

y1 := unapply(piecewise(0 <= t and t < t0, 0, y(theta(t))), t);

x2 := t->v2*t*cos(theta0);

y2 := t->v2*t*sin(theta0);

T := d*exp(theta0/sqrt(v1^2/v2^2-1))/(v1+v2);

A := animate(plot, [[x1(t), y1(t), t = 0 .. a], color = red, thickness = 4], a = 0 .. T, frames = 100);

B := animate(plot, [[x2(t), y2(t), t = 0 .. a], color = blue, thickness = 4], a = 0 .. T, frames = 100);

P1 := (x, y)->pointplot([[x, y]], color=red, symbol = solidcircle, symbolsize = 25);

P2 := (x, y)->pointplot([[x, y]], color = blue, symbol = solidcircle, symbolsize = 25);

A1 := animate(P1, [x1(t), y1(t)], t = 0 .. T, frames = 100);

B1 := animate(P2, [x2(t), y2(t)], t = 0 .. T, frames = 100);

display(B, B1, A, A1, scaling = constrained);

end proc:

Example:

Pursuit(1, 4, 1, 11*Pi*(1/6));

Workaround - use simplify with siderels option:

restart;

z := -d*(y+a)+a;

simplify(z, {a-a*d=c});

convert~(M,  parfrac);

For a more accurate picture of the graph of the function it is also useful to find the horizontal asymptote and the inflection points. The option scaling=constrained allows you to get the exact proportions of the plot:

f := x->(x^3-10*x^2-2*x+1)/(4*x^3+5*x+1):

VA := fsolve(4*x^3+5*x+1);   # The value x for vertical asymptote

HA := limit(f(x), x = infinity);   # The value y for horizontal asymptote

D1 := fsolve(numer((D(f))(x)) = 0);   # The zeroes of the first derivative

D2 := fsolve(numer(((D@@2)(f))(x)) = 0);   # The zeroes of the second derivative

Plots := plot([f(x), [VA, t, t = -3.3 .. 3.3], HA], x = -4.3 .. 4.3, -2.3 .. 2.3, discont, thickness = [2, 1, 1], color = [red, blue, blue], linestyle = [1, 3, 3]):   # The plots of the function and all its asymptotes

Points1 := plot(map(t->[t,f(t)], [D1]), style = point, color = green, symbol = solidcircle, symbolsize = 12): # Points extrema of the function (green color)

Points2 := plot(map(t->[t,f(t)], [D2]), style = point, color = yellow, symbol = solidcircle, symbolsize = 12):   # The inflection points of the function (yellow color)

plots[display](Plots, Points1, Points2, scaling = constrained);   # All in one plot

Variant with the point of tangency:

f := x->2*cos(x)-x^2: x0 := 2:

A := plot([f(x), f(x0)+(D(f))(x0)*(x-x0)], x = -1.3 .. 4.3, -17 .. 6, color = [red, blue], linestyle = [solid, dot]):

B := plot([[x0, f(x0)]], style = point, symbol = solidcircle, symbolsize = 17, color = blue):

plots[display](A, B);

Here is another solution, probably to be more clear:

restart:

with(LinearAlgebra):

A:=<3,0,0>: B:=<0,6,0>: C:=<0,0,4>: H:=<x,y,z>:

solve({DotProduct(B-A,C-H)=0, DotProduct(C-A,B-H)=0, Determinant(Matrix([A-H,B-H,C-H]))=0}):

'M'=convert(eval(H, %), list);

In the procedure  eul   in these lines

  x:= 0 + j*h:

  y := y[j]:

you assign new values to the formal parameters  x  and  y . It is invalid.