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.

@casperyc 

You can make it into column Vector by

V:= < v1, v2 >[..,1];

and the subs will work.

I don't think that it's strange that you want to make V a Vector. What's strange is that you want V1, ..., V4 to be Vectors also.

As to the "proper" way: I think it would be to define V first. Then if you really need V1, ..., V4 to be Vectors, create them from V like this:

s1:= a,b:
s2:= A,B:
s3:= seq(w[j],j=1..2):
s4:= w[3]*(x+y)^9:
V:= < s1, s2, s3, s4 >: or V:= Vector(7, [s1, s2, s3, s4]):
V1,V2,V3,V4:= V[1..2], V[3..4], V[5..6], V[7..7];

@casperyc 

You can make it into column Vector by

V:= < v1, v2 >[..,1];

and the subs will work.

I don't think that it's strange that you want to make V a Vector. What's strange is that you want V1, ..., V4 to be Vectors also.

As to the "proper" way: I think it would be to define V first. Then if you really need V1, ..., V4 to be Vectors, create them from V like this:

s1:= a,b:
s2:= A,B:
s3:= seq(w[j],j=1..2):
s4:= w[3]*(x+y)^9:
V:= < s1, s2, s3, s4 >: or V:= Vector(7, [s1, s2, s3, s4]):
V1,V2,V3,V4:= V[1..2], V[3..4], V[5..6], V[7..7];