Strange implementation (bug?) for numboccur.

restart;
L:=[1., HFloat(1.), 1., 0., 0.];
#                   L := [1., 1., 1., 0., 0.]

numboccur(L, 1.);
#                               2

numboccur(L, 1.0);
#                               0

numboccur(L, HFloat(1.0));
#                               1

evalb(L[1]=L[2]);
#                              true

select(`=`, L, 1.); select(`=`, L, 1.0);
#                          [1., 1., 1.]

#                          [1., 1., 1.]

dismantle(%);

#LIST(2)
#   EXPSEQ(4)
#      FLOAT(3): 1.
#         INTPOS(2): 1
#         INTPOS(2): 0
#      HFLOAT(2): 1.
#      FLOAT(3): 1.
#         INTPOS(2): 1
#         INTPOS(2): 0

arctand ignores the second argument.

arctand(y,x) simplifies to arctand(y)

BTW, for y <> 0,

arctan(y,x) = (Pi/2)*signum(y)-arctan(x/y)

The fact that leadterm is not documented under asympt is probably because asympt(f, x, n) is actually equivalent to series(f, x=infinity, n), and here leadterm appears. By taking n=2, leadterm is not strictly necessary, but it is convenient because O(...) is removed.

There is an old worksheet at the Application Center https://www.maplesoft.com/applications/view.aspx?SID=6322

It is clear enough and includes an animation.

Download extrema.mw

restart;
eq1:= f(x) + 2*f(1/x) + 3*f(x/(x-1)) = x:
eq:=eq1, simplify(eval(eq1, x=1/x)):
eq:=eq, simplify(eval(eq1, x=x/(x-1))):
eq:=eq, simplify(eval(eq1, x=-1/(x-1))):
eq:=eq, simplify(eval(eq1, x=1-x)):
eq:=eq, simplify(eval(eq1, x=(x-1)/x)):
solve([eq], indets([eq],function)): f(x) = eval(f(x), %);
simplify( eval(eq1,  f = unapply(rhs(%), x)) );   # check

So, you want the surface integral on the projection of a cuboid on a plane.

V:=<0,0,0>, <1,2,0>, <0,2,0>, <1,0,0>, <0,0,1>, <1,2,1>, <0,2,1>, <1,0,1>:
T:=<-1,0,0>, <0,-1,0>, <0,0,-1>: 
PM:=LinearAlgebra:-ProjectionMatrix([T[2]-T[1], T[3]-T[1]]):
Pr := v -> T[1]+PM.v:
Pr2:= v -> (PM.v)[[1,2]]:  # ConvexHull works in dimension 2.
P:=Pr~([V]):
P2:=Pr2~([V]):
ind:=ComputationalGeometry:-ConvexHull(Matrix(P2)^+):
C:=P[ind]: n:=nops(C):  # the vertices of the projection
Tri:=(v1,v2,v3) -> Surface( v1 + t*(v2+s*(v3-v2)-v1), s=0..1, t=0..1):
# the parametrization of the interior of a triangle
add(VectorCalculus:-SurfaceInt( x^2+y^2, [x,y,z] = Tri(C[1],C[i],C[i+1])), i=2..n-1);

Maple solved your equation with respect to z.
Note that LambertW is a function just like sin or exp.

Just increase Digits (exp(-200) is very small), e.g.

Digits := 20;

1. A faster method is to use Iterator.

with(Iterator):
C:=Combination(5,3):
K:=[1,2,3,4,a]:
seq( K[[seq(c[]+1)]], c=C);

        [1, 2, 3], [1, 2, 4], [1, 3, 4], [2, 3, 4], [1, 2, a], [1, 3, a], [2, 3, a], [1, 4, a], [2, 4, a], [3, 4, a]

2. To print the body of a procedure, use

interface(verboseproc=2);
print(combinat:-subsets);

or, 

showstat(combinat:-subsets);

evalf( Int(h, Omega, epsilon=1e-4, method=_CubaCuhre));

gives 4.12122491078804

(To have a probability you forgot probably to divide by (2*Pi)^3.)

You cannot write like this. diff(z, x)  simplifies to 0, because z is interpreted as an independent variable (just like x).
If z depends on x, you must write z(x). 
In this case,
diff(w(x,y,t,z(x)), x); 
==>  D[1](w)(x, y, t, z(x)) + D[4](w)(x, y, t, z(x))*diff(z(x), x) 

You may use e.g.
alias(z=z(x));
so that the argument x is omitted, but you must be careful (I personally do not use this facility in such situations).

You have a small mistake. Use:

Note that here the origin is at the minimal point of the catenary.