@Wolff Since the lower case l looks like an I in your title I just want to add that indeed the letter l has not been assigned a value. It appears in the denominator of Us.

For r<>0 there is no maximum for x=x(T) for T in the interval 0..500. 
If r > 0 then the limit of x(T) as T -> 0 (right) is infinity.

In plotting there is no need to use implicitplot as x is easily isolated (obviously):

restart;
Eq := x-(r+3.1*10^7*exp(-11616/(1.98*T)))/(3.1*10^7*exp(-11616/(1.98*T))+1.8*10^18*exp(-29691/(1.98*T))) = 0;
R := [0, 0.1e-2, 0.1e-1, .1, 1.6];

X:=solve(Eq,x);
plot([seq(X,r=R)],T=0..500,view=0..2);
limit(convert(X,rational),T=0,right);

@Carl Love The limit might do it, yes. Of course the assumption is that the expression is a linear combination of powers of t with coefficients independent of t.

@Carl Love Your power checks for type algebraic, but it won't work on a sum.

@Rouben Rostamian  My point was just that D doesn't produce a wrong result. It just returns unevaluated.

Just an observation. Nothing is broken in this sense:

convert(D(f)(x),diff);

So rather than a bug it might be called a weakness.

@Rouben Rostamian For your last plot the view found by plottools:-getdata is [1. .. 1., 0. .. 1.] .

It seems that this case of exact equality has been taken care of.

Continuing with the polar example:

evalindets(itv,float,round); #Now itv has 1 .. 1
plot(M,view=%);

@Christopher2222 A crude way would be to stop integration when the number of bounces has reach a certain level. That level, however, would be ad hoc. That is why I call it crude.
Here C__R = 0.8 and the level selected is 21.

sol := dsolve({deqs, ics,num(0)=0}, numeric, discrete_variables=[num(t)::integer],events = [[num(t)=21,halt],[y(t) = eval(surf, x = x(t)), [num(t)=num(t)+1,temp = diff(x(t), t), diff(x(t), t) = V_reflect[1], diff(y(t), t) = subs(diff(x(t), t) = temp, V_reflect[2])]]], output = listprocedure);

@darya In solving for the highest derivatives as in

SYS:=solve({sys},diff({A,C,F}(r),r,r) union {diff(B(r),r)});

I don't see any immediate problem.

Evaluating the right hand sides of SYS at r=ri using the initial conditions gives values that shouldn't frighten you. So I don't see what the problem is except that you surely have some huge numbers in the system.

@joshkalon tomleslie gave you the code.

I supply it here too.

restart;
ode:=diff(delta(x),x,x)+(1/2*0.3*(1+x)^3-0.7)/(0.3*(1+x)^3+0.7)/(1+x)*diff(delta(x),x)-3/2*0.3*(1+x)/(0.3*(1+x)^3+0.7)*delta(x)=0;
ics:=delta(0)=1,D(delta)(0)=-0.5;
dsolve({ode,ics});
res:=value(%) assuming x>0; #Presumably x>0 ?
eval(res,x=19); #Value at x=19
evalf(%); #In decimals
#Numerical check using numerical solution of the ode:
resnum:=dsolve({ode,ics},numeric);
plots:-odeplot(resnum,[x,delta(x)],0..19); #Plot of solution up till x=19
resnum(19); #Compares favorably with the result obtained from the exact solution.

Even for a flat earth you have the problem:

restart;
with(plots):
surf:=3;
V := Vector([diff(x(t), t), diff(y(t), t)]);
N := Vector([0,1]);
V_reflect := -(1+CR)*LinearAlgebra:-DotProduct(V,N,conjugate=false)*N+V;
CR := .8;
deqs := diff(y(t), t, t) = -9.81, diff(x(t), t, t) = 0;
ics := (D(x))(0) = 2, (D(y))(0) = 0, x(0) = 0, y(0) = 3.5;
sol := dsolve({deqs, ics}, {x(t), y(t)}, numeric,
   events = [[y(t) = 3, [temp = diff(x(t), t), diff(x(t), t) = V_reflect[1], diff(y(t), t) = subs(diff(x(t), t) = temp, V_reflect[2])]]], output = listprocedure);
xanim,yanim:=op(subs(sol, [x(t),y(t)]));
p1 := plot(surf, x = 0 .. 10, color = black):
animate(pointplot, [[xanim(t), yanim(t)], symbol = solidcircle, symbolsize = 15, color = black], t = 0 .. 3, frames = 150,background=p1);

By just using the initial condition pdsolve comes up with an answer, which happens also to satisfy the boundary condtitions:

pdsolve([pde,bc[3]]);

comes up with u(x, t) = x.

@fanicia Your list uses a capital X (upper case). In acer's reply it says x=X.

@J4James In this particular example you could do it in this (not very elegant) way:

res:=combine(PDEtools:-dchange({Xsubs2,Ysubs2,psubs2},Eq,[X,Y,P]));
indets(res,specfunc(anything,exp));
ln~(%) assuming real;
`=`(op(%));
simplify(%/epsilon);

@Markiyan Hirnyk The code exactly as I gave works for me in Maple 18, 17, 16, and 15.
So which Maple version are you using?

I don't have Maple 2015 yet.