There are several problems with your code. Firstly,  gamma  is a protected constant in Maple and you must either replace it with another symbol or execute  local gamma;  command. In addition, you have many parameters, so your symbolic expressions swell extremely fast and become completely useless. Therefore, I suggest that all (or almost all) parameters be given numeric values and calculate all numerically (by evalf command).

Example:

restart;
local gamma;
d(0):=1+1*I;
M(0):=1-2*I;
M__s:=1-I; d__s:=2+I; gamma:=3; nu:=4;
for k from 0 to 20 do 
u(k):=evalf(-d(k)/M(k));
p(k+1):=evalf(M__s*u(k)+d__s);
d(k+1):=evalf(d(k)+gamma*p(k+1)/(nu+abs(u(k))^2));
M(k+1):=evalf(M(k)+gamma*p(k+1)*conjugate(u(k))/(nu+abs(u(k))^2));
od;

I think that your problem can be effectively solved if the polynomial  P  is specified with numerical coefficients. To obtain a unique solution (of course, up to the sign of  Q  polynomial), we require that the degree of  R  does not exceed  n  (if the degree of  P  is 2*n) .

Example:

restart;
P := 9*x^4-6*x^3+13*x^2-8*x+5;
Q:=a*x^2+b*x+c:
R:=e*x+f:
T:=Q^2-R;
[seq(coeff(T,x^k)=coeff(P,x^k), k=1..4), tcoeff(T, x)=tcoeff(P)];
solve(%);
                  

PS. Of course, in this example, the polynomial  P is matched specifically to get integer coefficients. In general, as you can see, the problem reduces to solving a polynomial system of  2*n+1 equations with 2*n+1 unknowns.

This equation has only 1 positive root. Negative roots appear due to rounding errors. For negative values near 0, the function strictly decreases and tends to 0, remaining greater than 0 (this is easy to prove, but it is clearly seen from the plot):

restart;
fsolve(-2*sqrt(-K[1]^2+1)*K[1]+4*K[1]^2-Pi-2*arctan(sqrt(-K[1]^2+1)/K[1])=0, K[1]=0..1);
                                                0.9832314847
 

If I correctly understood your question, then the problem is as follows. There is a set of objects (in your example, these are matrices), each of which depends on a certain set of variables from the list var  (in your example these are  x  and  y ). Consider 2 objects and we raise  the question: whether one object be obtained from the other by some permutation from the list  var.

The procedure  IsEquivalent  solves this problem for matrices and any  var .

IsEquivalent:=proc(obj1::Matrix, obj2::Matrix, var::list)
local P, S;
uses combinat, LinearAlgebra;
P:=permute(var);
S:=[seq(subs(var=~p,obj2), p=P)];
if `or`(seq(Equal(obj1,s), s=S)) then true else false fi;
end proc:

Examples of use:

IsEquivalent(Matrix([[777,x,x],[77,7,7],[2,y,y]]), Matrix([[777,y,y],[77,7,7],[2,x,x]]), [x,y]);
IsEquivalent(Matrix([[777,x,x],[77,7,7],[2,y,y]]), Matrix([[777,y,y],[77,7,7],[2,x,y]]), [x,y]);

                                                               true
                                                              false

If you do work not with matrices, but with other objects (for example with polynomials), then a small obvious correction of this procedure will be required.

Never use square brackets to group summands. To do this, use only parentheses. The square brackets are used to create lists, and the braces ones to create sets.

Adjustment:

U:=eval(-(0.25*exp(1/2*epsilon*t)*(1/v^2-1+exp(epsilon*t))^(-0.5))*sin(4*t), [epsilon=1/10, v=2]);
plot(U, t=0..10, labels=[t,'U(t)'], size=[900,300]);

Solve separately for each function . B(t1,t2)  function can not be found from your system, so Maple returns NULL :

restart:
pde_sys :={diff(A(t1,t2),t1)*cos(B(t1,t2)) = 0, diff(A(t1,t2),t1)*sin(B(t1,t2)) = 0}:
dsolve(pde_sys, A(t1,t2));
dsolve(pde_sys, B(t1,t2));  # returns  NULL

                                         {A(t1, t2) = _F1(t2)}
                                                    

Edit.

 

The very existence of such a sequence is quite obvious without any formulas, because inequality  1 < (3^n) / (2^m) < 2  is equivalent to inequality  2^m < 3^n < 2^(m+1). In other words, any power of 3 is between some consecutive integer powers of 2  (which is absolutely obvious).

For any plotting, you must also specify the parameter values.

The example of a solution:

restart;
ode:=diff(u(t),t$2)-0.1*(1-64*u(t)^2)*diff(u(t),t$1)+16*u(t)=0;
Q1:=v/4*sin(4*t)+epsilon*(1/8*(v-v^3)*t*sin(4*t)+v^3/128*(cos(4*t)-cos(12*t)));
v:=1: epsilon:=1/2:
A:=DEtools:-DEplot(ode,u(t),t=0..10, [[u(0)=0, D(u)(0)=v]], u=-0.5..0.5):
C:=plottools:-getdata(A)[3]:
P:=plot(C, linestyle=dash, color=black, legend="Numerical Solution"):
B:=plot(Q1, t=0..10, u=-0.5..0.5,color=blue, legend="Regular Perturbation Expansion"):
plots:-display([P,B],title="Comparing the regular perturbation expansion to the numerical solution");

restart;
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:=unapply(Int(K1,z=0..1), F):
L1:=0.2:
d1:=0.2:
alpha:=Pi/6:
plots:-display(Vector[row]([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=[350,300]);

Output (vector-row of the plots with legends and titles):

Change the last line to following:

plots[display](seq(plots[display](q[i],r[i]), i=0..c-1), insequence=true, style=patchnogrid, labels=["x","y", "z"], scaling=constrained);
 

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.