Here is the solution with your data:

The code in 1d-math:

restart; 
rho := phi->chi/(1-e*cos(2*phi)):
e := 0.29:
chi := 0.5:
L := phi->Int(sqrt(rho(t)^2+(diff(rho(t), t))^2), t = 0 .. phi): 
plot(L(phi), phi = 0 .. 2*Pi); 
evalf(L(4)); 
evalf(L(5)); 
phi = fsolve(L(phi) = 50); 
plot(rho(phi), phi = 0 .. 2*Pi, coords = polar, scaling = constrained);
 

Download probléma_new.mw  # The code in 2d-math

Here are 2 more ways: 

# It does not even need in  simplify  command
eq1:=exp(I*phi);
conjugate(eq1)*eq1 assuming real;

# evalc command works under the assumption that all parameters are real
eq1:=exp(I*phi);
simplify(evalc(conjugate(eq1)*eq1));
 

restart;
interface(rtablesize=100):
for k from 1 to 10 do
A[k]:=<k,k^2; 2*k,2*k^2>;
od:
LinearAlgebra:-DiagonalMatrix([seq(A[k], k=1..10)]);

Probably you made syntax errors in  Rpm(Rpp(Psi00)) . For example, instead  e^(I*phi)  should be exp(I*phi) . e  is just a symbol in Maple.

Try the code below. And from now on, send a code in text form (which can be easily copied and then pasted into the Maple worksheet), not as a picture.

restart;
exp(-I*phi)*(-I*(2*I*exp(I*phi)*sqrt(1/Pi)*exp(-r^2/2)-2*I*exp(I*phi)*sqrt(1/Pi)*r^2*exp(-r^2/2))+2*exp(I*phi)*sqrt(1/Pi)*r*exp(-r^2/2)/r/(2*I*exp(I*phi)*sqrt(1/Pi)*r*exp(-r^2/2)))-2*r^2*exp(-I*phi)*exp(I*phi)*sqrt(1/Pi)*exp(-r^2/2);

simplify(%);
 

To find the global extrema in said region do:

f:=cos(x*y)*(x^2+y^2)^0.5:
minimize(f, x=0..1, y=0..1, location);
maximize(f, x=0..1, y=0..1, location);
plot3d(f, x=0..1, y=0..1);  # Visualization
 

Trunc := proc (eq::{polynom, procedure}, odr::nonnegint := 2, v::list := [x, y, z])
local a, b, q;
description " Truncates an algebraic equation to required degree";
a := eq; b := v;
map(select, q->evalb(degree(q, b) <= odr), `if`((eval(a))::procedure, a(b[]), a))
end proc:

A:=sqrt((1-b/r)^(-1));
B:=expand(op(1,A));
sqrt(``(expand(numer(B)*(-1)/``(denom(B)*(-1)))));
 

It seems that the syntax of  dsolve  command simply does not provide a vector input form. But the problem is easily solved by reducing to the standard syntax.

Example:

X:=t-><x1(t),x2(t)>:
A:=t-><1,3; 2,5>:
X0:=<1,2>:
Sys:=map(p->convert(p,list)[], {diff(X(t),t)=~A(t).X(t), X(0)=~X0});
dsolve(Sys, {x1(t),x2(t)});
 

r:=rand(-1...1.):

Examples of use:
X:=seq(r(), i=1..10);
Y:=seq(r(), i=1..10);

Addition:

In old versions of Maple do:

r:=()->RandomTools:-Generate(float('range' = -1 .. 1., 'method' = 'uniform')):

limit((u^m-v^m)/(u-v), u=v);
simplify(%);

Conversion to sum is easy to make for a specific  m , for example:
expand(factor((u^10-v^10)/(u-v)));

Addition. 1. Conversion to one side (from sum to closed form):
restart;
sum(u^(m-1-i)*v^i, i = 0 .. m-1);
simplify(%, symbolic);

2. For reverse conversion, you can use a special procedure (see the code and examples below). The expression must be entered in the form similar to  (u^m-v^m)/(u-v) . u, v  and  m   can be both symbols and/or numbers (if m is number then should be  m>1). To avoid premature calculations (in the case of numbers), pay attention to their correct input (see examples).

restart;
SpecialConv:=expr->Sum(expand(op([1,1,1],expr))^(op([1,1,2],expr)-1-i)*expand(op([1,2,2,1],expr))^i, i=0..op([1,1,2],expr)-1):

Examples of use:
SpecialConv((a^n-b^n)/(a-b));
SpecialConv((5^n-2^n)/``(5-2));
SpecialConv((a^10-b^10)/(a-b));
value(%);
SpecialConv((``(5)^10-``(2)^10)/``(5-2));
SpecialConv((a^10-``(2)^10)/(a-2));

Edit.

Try this:

RealInt := proc(f, x)
local A, B;
A:=int(f, x);
B:=applyrule(ln(a::anything)=ln(abs(a)), A);
simplify(B) assuming real;
end proc:

Examples of use:
RealInt(1+1/x, x);
RealInt(1/sin(x), x);
RealInt(x/(x^2+1), x);
RealInt(tan(x), x);
RealInt(diff(f(x),x)/f(x), x);
RealInt(1+1/abs(x), x);
RealInt(1/x, x) assuming x<0;
 

Eq1 := [(1/2)*(R[b]*kh[a2]+R[m]*kh[a2]-Rh[m]*kh[a2]+sqrt(R[b]^2*kh[a2]^2+2*R[b]*R[m]*kh[a2]^2-2*R[b]*Rh[m]*kh[a2]^2+R[m]^2*kh[a2]^2-2*R[m]*Rh[m]*kh[a2]^2+Rh[m]^2*kh[a2]^2))/kh[a2] = 0, 0 = 0, -(1/2)*(R[b]*kh[a2]+R[m]*kh[a2]-Rh[m]*kh[a2]+sqrt(R[b]^2*kh[a2]^2+2*R[b]*R[m]*kh[a2]^2-2*R[b]*Rh[m]*kh[a2]^2+R[m]^2*kh[a2]^2-2*R[m]*Rh[m]*kh[a2]^2+Rh[m]^2*kh[a2]^2))/kh[a2] = 0]:
simplify(Eq1) assuming kh[a2]>0, R[b]+R[m]-Rh[m]>=0 ;
simplify(Eq1) assuming kh[a2]<0, R[b]+R[m]-Rh[m]<=0 ;
                            [R[b] + R[m] - Rh[m] = 0, 0 = 0, Rh[m] - R[b] - R[m] = 0]
                            [R[b] + R[m] - Rh[m] = 0, 0 = 0, Rh[m] - R[b] - R[m] = 0]
 

We see that the expected result is correct if  kh[a2]  and  R[b]+R[m]-Rh[m]  have the same signs.
In the general case  in real domain  we have  sqrt(a^2*b^2)=abs(a*b) .

Edit.

restart;
f:=x->x^2:
A:=plot(f(x), x=0..2.5, color=blue, thickness=2, labels=[x, f(x)]):
T:=plots:-textplot([2.4,f(2), typeset("(",2,",",f(2),")")], font=[times,16]):
P:=plot([[2,f(2)]], style=point, color=red, symbol=solidcircle, symbolsize=15):
L:=plot([[2,t,t=0..f(2)], [t,f(2),t=0..2]], linestyle=2, color=black):
plots:-display(A, T, P, L, scaling=constrained, size=[300,500]);

The result:
                     

Edit.

If I correctly understood the question, then a simple procedure solves the problem:

restart;
f:=diff(ChebyshevT(n, r), r):
g:=unapply(simplify(eval(f, r = 0)), n):
seq(g(n), n=1..20);
piecewise(seq(op([n=k,g(k)]), k=1..20)); # In the form of piecewise

Addition. If, as indicated by Preben to use separate formulas for these cases (even  and  odd), then piecewise form can be written in a general form and shorter:

g:=n->piecewise(n::even, 0, n::odd, n*(-1)^((n-1)/2));
seq(g(n), n=1..20);  # Example of use

If we are interested in the set of all roots of some nonlinear equation, then the examples (OP's and vv's) show that we can not fully trust both  solve  command and fsolve and RootFinding:-Analytic  commands (the last two work in restricted intervals). Therefore, some qualitative analysis of the equation in question (even before its solution) is always useful, which concerns the existence and uniqueness of its roots. The simplest way to do this is to plot the corresponding function. Of course, the graph does not prove anything, but it allows you to notice some properties of the function, which you can then try to prove rigorously. In the example with OP's equation, it is easy to see the monotonicity of the function (it strictly increases). The analysis of the derivative shows that the function actually increases strictly on the whole real axis. Therefore, the equation can not have more than one real root. Since at the ends of the segment  [10, 20]  the function takes values ​​of different signs, the real root exists and it is unique:

Download root.mw