@Jalale  I did not understand what you mean by  "... animation with distances like AB/A'B' =AC/A'C'..."

@savalan I do not understand what you wrote at all. An equation does not depend on how you designate its roots. Your equations in the for loop are the same initial equations with the same roots, only the expression for the root  d  is shifted because you are actually making a change of variable in the same equations.

@savalan   What's the problem? Here are the first 50 positive values of  d:

seq(eval(rhs(sol[2]), _Z1=z), z=0..49);
0.03689590440, 0.5368959044, 1.036895904, 1.536895904, 2.036895904, 2.536895904, 3.036895904, 3.536895904, 4.036895904, 4.536895904, 5.036895904, 5.536895904, 6.036895904, 6.536895904, 7.036895904, 7.536895904, 8.036895904, 8.536895904, 9.036895904, 9.536895904, 10.03689590, 10.53689590, 11.03689590, 11.53689590, 12.03689590, 12.53689590, 13.03689590, 13.53689590, 14.03689590, 14.53689590, 15.03689590, 15.53689590, 16.03689590, 16.53689590, 17.03689590, 17.53689590, 18.03689590, 18.53689590, 19.03689590, 19.53689590, 20.03689590, 20.53689590, 21.03689590, 21.53689590, 22.03689590, 22.53689590, 23.03689590, 23.53689590, 24.03689590, 24.53689590

@quo  Yes of course! The functional operator  ->  in Maple is a special form of a procedure.

See help  ?Functional operators

@Markiyan Hirnyk 

1. I do not understand, what do you mean by "The parametrization L:=[[[t,sin(Pi*t/3)/2],t=0..3],[[3,0],[0,3],[1,1]], [[t,sqrt(t)],t=1..0]]: is done by hand". L is the actual parameter of the procedure and we always specify it by hand. Of course, in the text of the procedure it is easy to provide an explicit specification of some sections of the border instead of a parametric one. But I did not do this, so as not to complicate the text of the procedure by the different ways of specification of the border. Therefore, instead of [f(x), x=a..b]  we must write [[t, f(t)], t=a..b] .

2. The procedure returns an exact (symbolic) result. If Maple can not symbolically compute the integral, it returns it in an inert form and we can use numerical methods, for example  evalf  command. I think this approach is simpler than replacing the region of integration with a polygon.

See an answer to this question in  this post .

@acer You are a real virtuoso in such things!

Obviously, the original expression is periodic with the period  2*Pi , because It contains only the functions of  sin(phi)  and  cos(phi). Therefore, it is sufficient to solve the equation, for example, in the range  0 .. 2*Pi  and continue periodically to any desired interval. The same applies to finding extremums (local or global).

More precisely: combine(ee) shows that the period is  Pi .

Edit.

@Rouben Rostamian   Convert your comment in an answer and you will get a lot of votes up for it. Mine be the first.

For example you can write:

s := Int(exp(-x^2)*cos(2*x*y), x = 0 .. infinity):
s=value(s);

L:=[[0.5,0],[1,0],[1,7],[7,6.5],[5,9],[12,12],[10,13],[12,17],[10,13],[12,17],[8,16],[8,18],[3,16],[5,23],[3,21],[0.5,27]]:
L1:=ListTools:-Reverse(map(t->[1-t[1],t[2]],L)):
plot([L[ ],L1[ ]], color="#DC2828", thickness=5, axes=none, scaling=constrained);
# Or
plot([L[ ],L1[ ]], color="#DC2828", thickness=5, filled, axes=none, scaling=constrained);

@Franzs  Just now I noticed this new thread that duplicates your previous question. Look there my comment.

@Franzs   Maple shows the specified color only those planes that are present when setting the body. Add inequality  z>=-40 :

with(PolyhedralSets);
P1 := PolyhedralSet({-x-20 <= 0, -x-20 <= 20-y, -x-20 <= -(1/2)*y, -x-20 <= -z, -x-20 <= -(x+y)*(1/2), -x-20 <= -(y+z)*(1/2), z>=-40});
Plot(P1, axes = normal, view = [-40 .. 40, -40 .. 40, -40 .. 40], orientation = [30, -30, -60], color = blue, thickness = 2);
Plot(P1, axes = normal, view = [-40 .. 40, -40 .. 40, -40 .. -40], orientation = [0, 0, 0], color = blue, thickness = 2);

Rotation of the initial figure around the z axis:

restart;
with(PolyhedralSets):
P1 := PolyhedralSet({-x-20 <= 0, -x-20 <= 20-y, -x-20 <= -(1/2)*y, -x-20 <= -z, -x-20 <= -(x+y)*(1/2), -x-20 <= -(y+z)*(1/2)}):
A:=Plot(P1, axes = normal, view = [-40 .. 40, -40 .. 40, -40 .. 40], orientation = [30, -30, -60], color = blue, thickness = 2):
f:=alpha->plottools:-rotate(A, alpha, [[0,0,0],[0,0,1]]):
plots:-animate(plots:-display, [f(alpha)], alpha=0..2*Pi, frames=90);

@asa12   You wrote: If change to F(exp(t)) = 2*t, is it possible?

No, because for  t=exp(1)  we have:

F(F(exp(t))) = F(2*t) = F(2*exp(1)) = F(exp(ln(2)+1) = 2*(ln(2)+1) <> 0

@Markiyan Hirnyk   Use  expand  command first:

y=0.0000125698-0.0000125698*cos(54*x);
a:=convert(%,fraction);
b:=content(rhs(a), cos(54*x));
c:=primpart(rhs(a), cos(54*x));
y=b*``(c);
convert(%, float, 6);
diff(expand(rhs(%)), x);