restart; 
with(PDEtools):
umain := 1-exp(-y)+sum(A^i*u[i](y, z), i = 1 .. 5);
vmain := -1+sum(A^i*v[i](y, z), i = 1 .. 5);
wmain := sum(A^i*w[i](y, z), i = 1 .. 5);
pde[main] := diff(vmain, y)+diff(wmain, z);
pde[main[2]] := vmain*(diff(wmain, y))+wmain*(diff(wmain, z))-(diff(wmain, y$2));

sys[1] := [A^i*coeff(pde[main[2]], A) = 0, w[1](0, z) = sin(z), w[1](infinity, z) = 0];

for i from 2 to 5 do pde[i] := A^i*coeff(pde[main[2]], A^i) = 0; sys[i] := [pde[i], w[i](0, z) = 0, w[i](infinity, z) = 0] end do;

sol:=pdsolve(sys[1]);
w[1] := unapply(eval(w[1](y,z), sol), y,z);

# Examples of use
w[1](y, z);
w[1](1, 2);
plot3d(w[1], 0..1, 0..Pi);

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