Here are 2 ways of constructing the letter D. The first method is based on the standard method of constructing individual parts (circles, segments, ellipses), and then all together through  plots:-display  command. By this method we plot only the boundary, the interior remains unpainted.
The second method uses  Picture  procedure from here (for convenience, the code is presented here). This is a more automated way. As a parameter of the procedure, we indicate the list of individual parts of the border that are traversed by or counter-clockwise. In this way, we can paint the interior with the selected color and specify the boundary parameters.
The following are the codes for the standard worksheet. For a classic worksheet (for older versions Maple, there are some features for it), there is the link below.

restart;
Make_D:=proc(c::symbol:=black,h::posint:=2)
local B1, B2, C, E, F1, F2, F3, G1, G2, G3, H1, H2;
uses plots, plottools;
B1:=plot([-0.75+0.14*cos(t),0.49+0.14*sin(t), t=0..Pi/2], color=c, thickness=h):
C:=plot([0.79*cos(t),-0.13+0.8*sin(t), t=-Pi/2..Pi/2], color=c, thickness=h):
E:=plot([-0.14+0.66*cos(t),-0.13+0.717*sin(t), t=-1.94..1.94], color=c, thickness=h):
F1:=line([0,0.67],[-0.83,0.67], color=c, thickness=h):
F2:=line([-0.83,0.67],[-0.83,0.63], color=c, thickness=h):
F3:=line([-0.83,0.63],[-0.75,0.63], color=c, thickness=h):
H1:=line([-0.61,0.49],[-0.61,-0.76], color=c, thickness=h):
B2:=plot([-0.75+0.14*cos(t),-0.75+0.14*sin(t), t=-Pi/2..0], color=c, thickness=h):
H2:=line([-0.38,0.54],[-0.38,-0.8], color=c, thickness=h):
G1:=line([0,-0.93],[-0.83,-0.93], color=c, thickness=h):
G2:=line([-0.83,-0.93],[-0.83,-0.89], color=c, thickness=h):
G3:=line([-0.83,-0.89],[-0.73,-0.89], color=c, thickness=h):
display(B1, C,  E, F1, F2, F3, H1, B2, H2, G1, G2, G3, scaling=constrained);
end proc:

Examples of use:

plots[display](Make_D(red,3), plots[textplot]([1.6,0,"D"], font=[TIMES,310], color=grey), view=[-1..2.5,-1..1], size=[600,600]);  # For comparison, on the right, the same letter received by plots:-textplot  command
                

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:

Examples of use:

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(Picture(L2, color=white, [color=blue, thickness=5]), Picture(L1, color=yellow, [color=blue, thickness=5]), scaling = constrained, size=[400,400]):
LetterD;
                          

R:=alpha->plottools:-rotate(LetterD, alpha):
plots:-animate(plots:-display,['R'(alpha)], alpha=0..2*Pi, frames=60,  size=[600,600]);

LetterD1.mws

LetterD.mw
 

Edit.
 

I suggest using the functional style instead  eval : the code turns out to be more compact. This is also convenient if you want to calculate the values of  lambda  function at individual points. I also added text labels for your groups of graphs as in your picture:

h:=z->1-(delta2/2)*(1 + cos(2*(Pi/L1)*(z - d1 - L1))):
K1:=(4/h(z)^4)-(sin(alpha)/F)-h(z)^2+Nb*h(z)^4:
lambda:=(F,Nb,delta2)->Int(K1,z=0..1):

L1:=0.2:
d1:=0.2:
alpha:=Pi/6:
A:=plot( [seq(seq(lambda(F,Nb,delta2), Nb=[0.1,0.2,0.3]), F=[0.1,0.2,0.5])], delta2=0.02..0.1, color=[red$3,blue$3,black$3]):
B:=plots:-textplot([[0.05,-1,F=0.1],[0.05,1.5,F=0.2], [0.05,3,F=0.5]], font=[times, 14]):

plots:-display(A, B);

The output:

plot([4*t, 4*t-0.5*t^2, t=0..10]);

See  wiki  and the help on  ?plot,details

findOrderOf:=proc(g,p)
  local i:=1, gpwr:=1;
  for i from 1 to p-1 do
    gpwr:=gpwr*(g mod p);
    if gpwr = 1 then return i
    end if;
  end do;
end proc:

primes:=select(isprime,[`$`(2000 .. 3000)]);
select(p->findOrderOf(2,p)=p-1 and findOrderOf(3,p)=p-1,primes);

This works:

Download test_(1)_new.mw

In the body of the procedure  poly_out  replace the line

plots[display](trigr, lbase, trapezr, outAr, axes=none, scaling=constrained);

by the line

plots[display](trigr, lbase, trapezr, outAr, axes=normal, scaling=constrained);

We see that Maple reduces the calculation of this integral to a definite integral, which probably can not be calculated symbolically, but it is easily can be done numerically. Here is this calculation, as well as visualization: the cylindrical body is plotted, the volume of which is the computed integral:

Download Double_Integral.mw

restart;
q:=Matrix(3,10,(i,j)->i^2+j+2*i); 
with(plots):
V:=i->(q[..,round(i)]):
animate(arrow,['V'(k), width=[0.4,relative=false],head_width=[1,relative=false],color=red], k = 0.5..10.45, frames=60, axes=normal, paraminfo=false);
 

solve({7*cos(2*t)=7*combine(cos(t)^2)-5, t>=0, t<2*Pi}, t, allsolutions, explicit);  # Symbolic
evalf(%);  # Numeric
                 

Let the function  f  be the inverse of the function temperature vs chemicals, M  is your Array. Do the following:

A:=convert(M, Matrix);
L:=convert(<f~(A[..,2]) | A[..,1]>, listlist);
plot(L);
 

eq16:=r(t)=d[vol]*V/(KUS*V^2+L*tau);
r(t)=(numer(rhs(eq16))/V)/'(expand(denom(rhs(eq16))/V))':
eq17:=%;

Slightly reducing the number of digits, it is easy to check the identity of the lists:

evalb(evalf[8](listA) = evalf[8](listB));
 

restart;
f:=unapply(<r(t)*cos(t), r(t)*sin(t)>, t):
v:=diff~(f(t), t);  # The vector of velocity

Example of use. The calculation of the vector of velocity for the curve  r(t)=t  (Archimedes' spiral) at the point  t=1 :

f:=unapply(<t*cos(t),t*sin(t)>, t):
v:=unapply(diff~(f(t),t), t):
v(t);  # The vector of velocity
v(1);  # Velocity at the point  t=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), scaling=constrained);  # Visualization

Edit.

Obviously, this is a bug. Here's a very trivial example:

f1 := x->x^2:   f2 := x->x^2:
is(f1=f2);
                                                false

Probably  is  command simply does not work with procedures. Obviously to check the equality of the two procedures  f  and  g , it suffices to verify that  f(x)=g(x)  for an arbitrary  x:

f1:=x->x^2:  f2:=t->t^2:
is(f1(a) = f2(a));
                                                 true

Use  evalc  command for this.

Example:

evalc(exp(k*Pi*I));
                                       cos(k*Pi)+I*sin(k*Pi)