restart;
sol := dsolve({(diff(r*(diff(u(r), r)), r))/r = 0, u(R0) = U2, u(h) = U1}, u(r));
plot(eval(rhs(sol), [h = 1, R0 = 5, U1 = 2, U2 = 3]), r = 1 .. 5);

It is much easier to calculate the length of this arc of the ellipse (half of the total ellipse) using its parametric equations:

eqn := (1/32)*(x-16)^2+(1/2025)*(y+0)^2 = 1;
a:=sqrt(32):  b:=sqrt(2025):
int(sqrt(diff(16+a*cos(t),t)^2+diff(b*sin(t),t)^2), t=0..Pi);
evalf(%);
plots[implicitplot]([(1/32)*(x-16)^2+(1/2025)*(y+0)^2 = 1], x = 0 .. 32, y = 0 .. 60, scaling = constrained);

To plot arrows, use  plots:-arrow  command, and for text use  plots:-textplot command. See help on these commans for details.

Example of the plotting:

with(plots):
P:=plot([seq(5-a/20*(x+5)^2, a=[1,2,3])], x=-5..3, -2..7, color=[red,blue,green], size=[700,300]):
A1:=arrow([-2,6],[-0.5,-1.85], width=[0.01,relative=false], head_width=[0.07,relative=false], head_length=[0.2,relative=false]):
A2:=arrow([-1.3,6],[-0.5,-1.9], width=[0.01,relative=false], head_width=[0.07,relative=false], head_length=[0.2,relative=false]):
A3:=arrow([-0.6,6],[-0.5,-1.65], width=[0.01,relative=false], head_width=[0.07,relative=false], head_length=[0.2,relative=false]):
T:=textplot([[-0.6,6.4,a=1],[-1.3,6.4,a=2],[-2,6.4,a=3]], font=[times,14]):
display(P,A1,A2,A3,T);

Your integral contains a large number of parameters. Integral is easy to calculate if you specify the values of these parameters. Also note that you should not use a variable and an indexed variable  with the same name in the same worksheet (you have d and d[i], x and x[i]).

Example:

restart;
eval(sum(d[i]*ln(x[i]/theta), i = 1 .. 3)/(d1*(sum(ln(x[i]/theta), i = 1 .. 3))^v),[v=2,d1=6,d[1]=1,d[2]=2,d[3]=3,x[1]=1.5,x[2]=2.3,x[3]=3.7]);
int(%, theta=0..2);
                          

Try this version. In classic GUI of Maple 2017.3 it works:

restart;
with(plots):
with(plottools):
w:=1.5:
lup := line([w,0.9], [w,5.1], color=red, thickness=4):
plt:=plots[implicitplot]((x-1)^4/24+(y-3)^4/10=2, x=w..7, y=0..12,color=red, thickness=4):
op([1,1],plt);
convert(%,listlist):
P:=polygon(%, color=gold):
plots[display]([lup,plt],scaling=constrained);
plots[display]([lup,P,plt],color=gold,scaling=constrained);

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);