Using evalf is risky: If the precision is not sufficiently high, you may get a wrong result---not just the absence of a result. Whereas is will add digits of precision sufficient to determine the inequality with certainty. For example:

Digits:= 3:
if Pi < 22/7 then Yes else No fi;
Error, cannot determine if this expression is true or false: Pi < 22/7

if is(Pi < 22/7) then Yes else No fi;
                              Yes
evalf(Pi), evalf(22/7);
                           3.14, 3.14
if evalf(Pi) < evalf(22/7) then Yes else No fi;
                               No

While this example may seem contrived, it is very common in numerical programming for an inequality to be decided by the last digit of precision.

Using evalf is risky: If the precision is not sufficiently high, you may get a wrong result---not just the absence of a result. Whereas is will add digits of precision sufficient to determine the inequality with certainty. For example:

Digits:= 3:
if Pi < 22/7 then Yes else No fi;
Error, cannot determine if this expression is true or false: Pi < 22/7

if is(Pi < 22/7) then Yes else No fi;
                              Yes
evalf(Pi), evalf(22/7);
                           3.14, 3.14
if evalf(Pi) < evalf(22/7) then Yes else No fi;
                               No

While this example may seem contrived, it is very common in numerical programming for an inequality to be decided by the last digit of precision.

My first answer did not address the issue of coloring the edges. That is a little bit more complicated. Here are the faces and edges together.

restart:
rand256:= rand(0..255):
RandColor:= ()-> COLOR(RGB, evalf(['rand256'()/256 $ 3])[]):
C:= plottools:-cuboid([0$3], [1$3], color= ['RandColor'() $ 6]):
faceplot:= plots:-display(C, transparency= 0.01, style= patchnogrid, glossiness= 1):

Do some combinatorics to generate the edges. First generate the vertices.
V:= ['('([i,j,k] $ i= 0..1)' $ j= 0..1)' $ k= 0..1];
    [[0, 0, 0], [1, 0, 0], [0, 1, 0], [1, 1, 0], [0, 0, 1],  [1, 0, 1], [0, 1, 1], [1, 1, 1]]

Now the edges: From all pairs of points, select those unit distance apart.
E:= select(e-> LinearAlgebra:-Norm(< `-`(e[]) >, 2) = 1, combinat:-choose(V,2));
     [[[0, 0, 0], [0, 0, 1]], [[0, 0, 0], [0, 1, 0]],  [[0, 0, 0], [1, 0, 0]], [[0, 0, 1], [0, 1, 1]],
     [[0, 0, 1], [1, 0, 1]], [[0, 1, 0], [0, 1, 1]],  [[0, 1, 0], [1, 1, 0]], [[0, 1, 1], [1, 1, 1]],
     [[1, 0, 0], [1, 0, 1]], [[1, 0, 0], [1, 1, 0]],  [[1, 0, 1], [1, 1, 1]], [[1, 1, 0], [1, 1, 1]]]

edgeplot:= plots:-display(
     [seq](plottools:-line(e[], color= RandColor()), e in E),
     thickness= 5, transparency= 0.45, glossiness= 0
):
plots:-display([edgeplot, faceplot], lightmodel= light4, axes= none);

My first answer did not address the issue of coloring the edges. That is a little bit more complicated. Here are the faces and edges together.

restart:
rand256:= rand(0..255):
RandColor:= ()-> COLOR(RGB, evalf(['rand256'()/256 $ 3])[]):
C:= plottools:-cuboid([0$3], [1$3], color= ['RandColor'() $ 6]):
faceplot:= plots:-display(C, transparency= 0.01, style= patchnogrid, glossiness= 1):

Do some combinatorics to generate the edges. First generate the vertices.
V:= ['('([i,j,k] $ i= 0..1)' $ j= 0..1)' $ k= 0..1];
    [[0, 0, 0], [1, 0, 0], [0, 1, 0], [1, 1, 0], [0, 0, 1],  [1, 0, 1], [0, 1, 1], [1, 1, 1]]

Now the edges: From all pairs of points, select those unit distance apart.
E:= select(e-> LinearAlgebra:-Norm(< `-`(e[]) >, 2) = 1, combinat:-choose(V,2));
     [[[0, 0, 0], [0, 0, 1]], [[0, 0, 0], [0, 1, 0]],  [[0, 0, 0], [1, 0, 0]], [[0, 0, 1], [0, 1, 1]],
     [[0, 0, 1], [1, 0, 1]], [[0, 1, 0], [0, 1, 1]],  [[0, 1, 0], [1, 1, 0]], [[0, 1, 1], [1, 1, 1]],
     [[1, 0, 0], [1, 0, 1]], [[1, 0, 0], [1, 1, 0]],  [[1, 0, 1], [1, 1, 1]], [[1, 1, 0], [1, 1, 1]]]

edgeplot:= plots:-display(
     [seq](plottools:-line(e[], color= RandColor()), e in E),
     thickness= 5, transparency= 0.45, glossiness= 0
):
plots:-display([edgeplot, faceplot], lightmodel= light4, axes= none);

@somayeh You're welcome. Let me know if you are having any trouble using my answer.

@somayeh You're welcome. Let me know if you are having any trouble using my answer.

Yes, I should have added that my identity held for positive a and b.

Is there some reference where one can find all codes of the form `&...;`? Is this part of some system of codes that is more universal than Maple?

The following works in 1D Input, typed just as is, no mouse work required (the way I like it):

plot(
     [x, 2*x], x= -1..3,
     axis[1]= [gridlines= [[1], linestyle= dot]],
     tickmarks= [[1=a, 0= "0"], [1=b[z]^`&lowast;`, 0=0]],
     legend= ["", b[z]^`&lowast;`],
     axesfont= [TIMES,ROMAN,16],
     legendstyle= [font= [TIMES,ROMAN,16]]
);

The difference is that b_z becomes b[z].

Is there some reference where one can find all codes of the form `&...;`? Is this part of some system of codes that is more universal than Maple?

The following works in 1D Input, typed just as is, no mouse work required (the way I like it):

plot(
     [x, 2*x], x= -1..3,
     axis[1]= [gridlines= [[1], linestyle= dot]],
     tickmarks= [[1=a, 0= "0"], [1=b[z]^`&lowast;`, 0=0]],
     legend= ["", b[z]^`&lowast;`],
     axesfont= [TIMES,ROMAN,16],
     legendstyle= [font= [TIMES,ROMAN,16]]
);

The difference is that b_z becomes b[z].

Complex arithmetic allows for this transformation:

ln(a) - ln(b) = ln(a/b) = ln(-a/-b) = ln(-a) - ln(-b)

Now consider your eq16:

Arranging the ln terms on the left side and the non-ln terms on the right and factoring out hmax, we get on the left side

hmax*(ln(-hmax) - ln(h - hmax))

Applying what I said in the first sentence, this is equivalent to

hmax*(ln(hmax) - ln(hmax - h)),

where now all the arguments of ln are positive.

This is reminiscent of ?evalr and shake----range arithmetic.

A problem still arises if you want a local unevaluated name gamma. There is nothing that I can do to prevent evalf(gamma) from returning a number. Even `evalf/gamma`:= NULL doesn't seem to help.

A problem still arises if you want a local unevaluated name gamma. There is nothing that I can do to prevent evalf(gamma) from returning a number. Even `evalf/gamma`:= NULL doesn't seem to help.

@Preben Alsholm `?()` is the function invocation operator. Think of the expression f(x) as being a binary infix operation whose operands are f and x. If you want to treat that operator as a name, it is `?()`. The most documentation that I can find about it is at ?use , at the end of "Description" where the overloadable operators are described.

`?()`(f, [x,y]) is the same as f(x,y). I had not thought of apply. Yes, the are almost the same, the difference being that `?()` requires the arguments in square brackets. Also, `?()` is kernel and apply is library.

@Preben Alsholm `?()` is the function invocation operator. Think of the expression f(x) as being a binary infix operation whose operands are f and x. If you want to treat that operator as a name, it is `?()`. The most documentation that I can find about it is at ?use , at the end of "Description" where the overloadable operators are described.

`?()`(f, [x,y]) is the same as f(x,y). I had not thought of apply. Yes, the are almost the same, the difference being that `?()` requires the arguments in square brackets. Also, `?()` is kernel and apply is library.