Similar problems arise only in the classic interface. In standard one polygon's filling works correctly for any simply-connected polygons (convex or concave). I added only one line to your code (in Maple 2017.3):

plots[display](<l | lM | H>,scaling=constrained, color=gray);

@das1404  The file is so big because of the animation. Try not to open this file, but just copy the following commands as text and run them in your Maple. I think this should work in Maple 7 (there are no references to the font here). Run not all the commands at once, but one at a time.

restart;
Picture := proc (L, C, N::posint := 100, Boundary::list := [linestyle = 1]) 
local i, var, var1, var2, e, e1, e2, P, Q, h; 
global Border; 
for i to nops(L) do 
if type(L[i], listlist(algebraic)) then P[i] := op(L[i]) else 
var := lhs(L[i, 2]); var1 := lhs(rhs(L[i, 2])); var2 := rhs(rhs(L[i, 2])); 
h := (var2-var1)/N; 
if type(L[i, 1], algebraic) then e := L[i, 1]; 
if nops(L[i]) = 3 then P[i] := seq(subs(var = var1+h*i, [e*cos(var), e*sin(var)]), i = 0 .. N) else P[i] := seq([var1+h*i, subs(var = var1+h*i, e)], i = 0 .. N) end if else 
e1 := L[i, 1, 1]; e2 := L[i, 1, 2]; P[i] := seq(subs(var = var1+h*i, [e1, e2]), i = 0 .. N) end if end if end do; 
Q := [seq(P[i], i = 1 .. nops(L))]; 
Border := plottools[curve]([op(Q), Q[1]], op(Boundary)); 
[plottools[polygon](Q, C), Border]; 
end proc:

L1:=[[[0.79*cos(t),-0.13+0.8*sin(t)], t=-Pi/2..Pi/2], [[0,0.67],[-0.83,0.67],[-0.83,0.63],[-0.75,0.63]], [[-0.75+0.14*cos(t),0.49+0.14*sin(t)], t=Pi/2..0], [[-0.61,0.49],[-0.61,-0.76]], [[-0.75+0.14*cos(t),-0.75+0.14*sin(t)], t=0..-Pi/2], [[-0.75,-0.89],[-0.83,-0.89],[-0.83,-0.93],[0,-0.93]]]:
L2:=[[[-0.14+0.66*cos(t),-0.13+0.717*sin(t)], t=-1.94..1.94], [[-0.38,0.54],[-0.38,-0.8]]]:

LetterD:=plots[display](plot([0.62,[-0.75+0.14*cos(t),0.48+0.14*sin(t),t=0..Pi/2]], x=-0.83..-0.75, color=white,filled), Picture(L2, color=white, [color=blue, thickness=5]), Picture(L1, color=yellow, [color=blue, thickness=5]), scaling = constrained):
LetterD;

N:=60: # The number of the frames
R:=seq(plottools[rotate](LetterD, 2*Pi*k/N), k=0..N):
plots[display](R, insequence=true);

@das1404  You can see here an example of creating and rotating in 2D the letter D. It seems you did not see this my answer. I'm wondering whether this file   LetterD1.mws  will work in Maple 7?

@waseem   Change the last line to

plots:-display(Vector([seq(plot([seq(eval(lambda(F), Nb=j), j=[0.1,0.2,0.3])], delta2=0.02..0.1, legend=[Nb=0.1,Nb=0.2,Nb=0.3], title=('F'=F)), F=[10,20,30])]), size=[900,500]);

@ganelon  From this specification nothing will change: the solution is unique for any  y>0

@acer  How did you know about [Gr,Br] ?

@Carl Love  You're right.

@das1404   If I need to draw several graphs on one plot, then I specify the common options (axes, scaling, view) only at the end of  plots[display]  command. In this case there will be no contradiction between the options.

@Preben Alsholm 

We can simplify a little more by

EXPR2:=combine(EXPR1, ln) assuming h>0, R0>0;
length(EXPR2);

@alfarunner  There are no problems here. The coefficients before rho hat and phi hat are the coordinates of the vector v(t) in this basis (here phi=t):

restart;
f:=unapply(<r(t)*cos(t),r(t)*sin(t)>, t):
v:=unapply(diff~(f(t),t), t):
V:=<cos(t),sin(t); -sin(t),cos(t)>.v(t);
simplify(V);  # The coordinates in the new basis (rho hat and phi hat)

I have never worked in the Physics package.

@alfarunner  Above I wrote how to calculate the coordinates of the velocity vector in standard Cartesian coordinates. If you need the same, but in other coordinates, you just need to do a coordinate transformation, multiplying by the corresponding matrix:

restart;
f:=unapply(<r(t)*cos(t),r(t)*sin(t)>, t):
v:=unapply(diff~(f(t),t), t):
V:=<cos(t),sin(t); -sin(t),cos(t)>.v(t);

f:=unapply(<t*cos(t),t*sin(t)>, t):  # The same example
v:=unapply(diff~(f(t),t), t):
v(t);
v(1);
simplify(<cos(1),sin(1); -sin(1),cos(1)>.v(1));
plots:-display(plot([convert(f(t),list)[],t=0..3], color=red, thickness=3), plots:-arrow(f(1),v(1), width=[0.01,relative=false],head_width=[0.08,relative=false],color=blue), plots:-arrow(f(1),<cos(1),sin(1)>, width=[0.01,relative=false],head_width=[0.08,relative=false],color=black), plots:-arrow(f(1),<-sin(1),cos(1)>, width=[0.01,relative=false],head_width=[0.08,relative=false],color=black), scaling=constrained);
 

@Preben Alsholm  Thank you for clarifying. Yet from the standpoint of pure mathematics it seems that this is not a very good design.

@minhthien2016  To do this, the code should be rewritten as a procedure. After that, you can apply this procedure to any cube by specifying as a procedure argument the list of cube vertices.

Here is this procedure:

restart;
RegTetr:=proc(L)
local T, Dist, P;
uses combinat;
T:=combinat:-choose(L, 4);
Dist:=(p,q)->sqrt((p[1]-q[1])^2+(p[2]-q[2])^2+(p[3]-q[3])^2);
P:=v->combinat:-choose(v, 2);
[seq(`if`(nops(convert((Dist@op)~(P(t)),set))=1, t, NULL), t=T)];  
end proc:

Example of use:

L:=[[0,0,0], [1,0,0], [1,1,0], [0,1,0], [0,0,1], [1,0,1], [1,1,1], [0,1,1]]:
L1:=2*~(map(t->t-~[1/2,1/2,1/2], L));
RegTetr(L1);

@vv  Your solution is impressive, but weakly motivated (much remains behind the scenes). Why did you take these ranges  x=50..100, y=30..100, H=400..500 ? In addition, we can not so round up  y=40.00000000, because we obtain 0 under the logarithm sign.

Here is another solution in which we first eliminate  H , simplify the obtained system , and then decide, taking as the initial point x=41, y=41, because obviously should be  x>40, y>40  for real solution :

restart;
Digits:=1000:
Sys:=[21000 = ((10/1000)/60)*4181*983*(x-y),
H = -((10/1000)/60)*4181*ln(1-((x-y)/(x-40))),
H = 124+102*ln(x-40)+(7+2*ln(x-40))];
A:=eliminate(Sys, H):
Sys1:=convert(A[2], list):
Sys2:=expand(Sys1) assuming x>40, y>40;
fsolve(Sys2, {x=41,y=41}):
Sol:={%[], eval(Sys[3],%)};  # The final result
eval(Sys, Sol):  # Check
evalf[20](%);
 

Download System.mw

@Aaeru Michi  Maple 16 just does not cope with this example (I solved in Maple 2017). Here are 2 ways that maybe help you:

minimize(V, x = 0 .. 1, location) assuming x>0;  # or

minimize(-(1/2)*abs(x)*(c*x^2-2*L+x)/(2*c*x+1)^2, x = 0 .. 1, location);  # I just removed  signum(2*c*x+1)  from your example