@Markiyan Hirnyk  Thank you.

I can not agree with your statement  "This approach is rather a handicraft than a mechanization because a parametrization of the boundary of a domain  is needed."  And if some region is initially defined by the equation of its boundary? In addition, not every region can be written as the intersection of individual regions. For example, a nonconvex polygon can not be defined in this way, because the intersection of half-planes is always convex.

@maxburakho  See Help on  ?fsolve,detail . I did not find there any options that would allow me to force a solution method. In fact, fsolve command itself chooses the most suitable method for calculating the root. For your case, the system is actually solved exactly (not by Newton method), solve command gives the same result. Probably you will have to write code for solving the system yourself, using Newton's method or to look for something similar on the Internet.

@lucaud  I did not understand your request regarding volume. The volume is already calculated in Carl's answer (2 ways) and in my answer (1 way).

When creating complex animations, it is more convenient and quicker to use a procedure to create one frame for an arbitrary value of the animation parameter, and then output successive frames with plots:-animate command.

Here is this technique for your example. I also increased the number of frames for better perception and change the color, otherwise it looks rather gloomy:

restart;
Australopithecus := [[75, 25], [97, 30], [93, 40], [93, 45], [83, 50], [80, 55], [79, 60], [81, 73], [74, 76], [68, 81], [60, 82], [50, 83], [40, 80], [30, 71], [25, 60], [24, 50], [25, 37], [15, 33], [10, 30], [45, 10], [55, 16], [65, 10], [80, 8], [93, 14], [96, 24]]:
man := [[95, 39], [113, 40], [111, 47], [118, 53], [113, 62], [109, 72], [112, 88], [112, 95], [107, 112], [99, 117], [85, 122], [72, 122], [49, 117], [36, 104], [31, 78], [39, 52], [43, 43], [44, 34], [39, 16], [73, 3], [81, 17], [98, 14], [105, 17], [104, 26], [111, 33]]:

morph := proc(poly1, poly2, t)
uses plots, plottools; 
if nops(poly1) <> nops(poly2) then ERROR("mensaje.") end if; 
display(polygon((1-t)*poly1+t*poly2, color = gold), axes = none); 
end proc:

plots:-animate(morph, [Australopithecus, man, t], t = 0 .. 1, frames = 60, paraminfo = false);

Formulate more exactly what you want: volume, picture or something else. Specify these 2 solids here.

Here is a procedure based on vv's idea:

restart;
LR:=proc(X::list)
local X1, y0, X0, L, x0, Y0, y1;
X1:=simplex:-convexhull(X);
y0:=min(seq(x[2], x=X1));
X0:=sort(select(t->t[2]=y0, X1),key = (x->x[1]));
L:=select(t->t[1]>=X0[-1,1], X1);
x0:=max(seq(l[1], l=L));
Y0:=sort(select(t->t[1]=x0, L),key = (x->x[2]));
y1:=min(seq(y[2], y=Y0));
remove(t->t[2]>y1, L);
end proc:

Example of use:
LR([[1,1],[1,2],[1,3],[2,1],[3,2],[3,6],[4,2],[4,6],[5,4],[5,5]]);
                  [[2, 1], [4, 2], [5, 4]]

Visualization of this example:
A:=plot([[1,1],[1,2],[1,3],[3,2],[3,6],[4,6],[5,5]], style=point, color=red, symbol=circle, symbolsize=20):
B:=plot([[2,1],[4,2],[5,4]], style=point, color=red, symbol=solidcircle, symbolsize=20):
P1:=plot([[2,1],[4,2],[5,4]], color=blue, thickness=2):
P2:=plot([[5,4],[5,5],[4,6],[3,6],[1,3],[1,1],[2,1]], linestyle=3, color=blue):
plots:-display(A,B,P1,P2, view=[0..5,0..6], scaling=constrained);

The code was edited.

@Earl  As a point corresponding to the initial value of the parameter t=0, we can take any point of this circle. As you can see from the text of the code, first I found the center of this circle (the list C=[1,1/2,1/2]), so the plane passing through the center intersects the circle at least in two points. Solving the corresponding system, I found 2 points and as a point corresponding to t = 0, I took a point with a larger first coordinate.

@max125  It's like a staircase:

plot3d(trunc(x), x=0.99..4.99,y=0..5);

@Markiyan Hirnyk I do not understand what is the point of getting a formal analytical expression without  checking its correctness and how all this will work in special cases?

Please see:

restart; 
Expr:=value(int((x-a)^m/x, x = a .. x)) assuming x > a, a > 0, m>-1;
eval(Expr, [a=2,m=3]);
eval(Expr, [a=-2,m=3]);
eval(Expr, [a=2,m=-3]);
 

We get continuous errors. Can you get, using  Expr, the same results as in my answer above?

@vv  Thank you! Corrected. When I was adding comments to the code, I accidentally deleted one line of code. Now everything works properly.

How did you find this value  ε=155437 / 1250000? In fact, it is overestimated (about 30 times). You can find it in Maple as follows:

f:= x*(1 - cos(x))/(x - sin(x)):
L:=limit(f, x=0);
epsilon:=maximize(abs(f - L), x=-0.2..0.2);
(155437/1250000)/epsilon;

                                    3
                         0.004000393186
                           31.08434452
 

@brian bovril  No, I do not have such a program. In fact, the reprint requires no more than 5-7 minutes, because the first 3 equations are very similar. I typed one equation, and then copied and pasted it 2 times, and then made the necessary changes.

@Markiyan Hirnyk  

1. I understand it. Of course, it would be better if Maple imposed the restriction  _Z1>=0  and  _Z1<=n-1

2. What is the point of asking Maple to write a formula that can be found in any textbook. Maple is good in specific calculations, which was shown in my answer.

@Markiyan Hirnyk  If we follow your (general) solution, then we get (formally) infinitely many roots:

Originally _Z1, renamed _Z1~:
  is assumed to be: integer

@max125  You missed the semicolon (or the colon) after
  a:=solve(Sys, [x], allsolutions, explicit=true)