solve(3^665/2^x=y, x);

For a more compact output (to prevent exponentiation  3^665) we can do so:

solve(3^``(665)/2^x=y, x);
                                        

Edit.

We put in correspondence to each point  (x,t)  the vector  <Re(u(x,t),), Im(u(x,t))> . To plot the corresponding vector field, we can use  plots:-fieldplot  command. The length of each variable vector is the modulus of the corresponding complex number:

Download plot33_b.mw

Never use in one formula an indexed variable , say  r[i] , and just a variable  with the same name  r, because they are logically connected. For example, if  r  is a list, then  r[i]  is  i_th  element of this list. 

As an integration variable for a definite integral, take a symbol that is not in your document, for example  _s :

int(_s*r[bc]*r[tc], _s = r[bc] .. r[tc]);

RootOf  is a method of compact implicit representation of the roots of an equation. To explicitly represent the roots, allvalues  command is used (see the example below).

solve(x^4-3*x-1=0);
allvalues([%]);

The last line of your code (in  b ) does not already have RootOf , all the roots are in explicit form.

See help on  ?RootOf  for details.

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:=%;