@Carl Love When trying your code in Maple 2015 exactly as written I got

Error, cannot determine if this expression is true or false: FAIL < 10000

I guess you left out setting Digits:=15 or similar, as I just noticed that the code works then.

Have a look at the help page
?pdsolve, numeric

In particular take a look at the last example at the bottom of the page.

@JohnPo To get data from the plots you can use plottools:-getdata (see below). However the module returned by pdsolve/numeric exports the procedure value, which may be what you want.

pds := pdsolve(eval(PDE,s=1), IBC, numeric);
## Using value:
val:=pds:-value(theta(y,z)); #theta(y,z) can be left out
##val is a numerical procedure taking two arguments, y and z:
val(0.1234,2.1);
# Alternatively use option listprocedure:
pds:-value(theta(y,z),output=listprocedure);
T:=subs(%,theta(y,z));
T(0.1234,2.1);
########## Getting data from plots:
pds:-plot3d(z=0..3); p1:=%:
plottools:-getdata(p1);
A:=%[3];
## A[i,j] contains the theta value at the (y,z) point corresponding to the grid point (i,j) in the 25x25 equally spaced grid over the yz-rectangle [0,1]x[0,3]:
plots:-matrixplot(A);
#1D plot:
pds:-plot(z=3); p2:=%:
plottools:-getdata(p2);
B:=%[3];
#The meaning of B is different from the 2D-case: First column contains the y's, second the theta's.
B[1,..];
plot(B);

@JohnPo It is rather doubtful with the initial and boundary conditions that pdsolve can give you a symbolic solution.

Try without the conditions. Begin by removing the list brackets around PDE for convenience.

PDE := (1-y^(s+1))*(diff(theta(y, z), z)) = diff(theta(y, z), `$`(y, 2))+y^(s+1)*(s+1)^(1/s+1);
pdsolve(PDE); #You see separation of variables being used
pdsolve(PDE,build); #The odes somewhat solved: Notice the DESol structure
# Now a direct attempt at including IBC:
#Another syntax for inputting IBC is used in the symbolic case:
pdsolve({PDE} union IBC); #IBC was defined as a set.

Failure, however

@acer As one who always uses the Standard GUI with 1D-input I was surprised to learn that the necessity for colon or semicolon at the end of a statement had been removed.
I don't see the need for that at all.
Has it been done to help relatively inexperienced Maple users? Wouldn't they in fact be using 2D-input? (Not that I would recommend that, though).

@ramakrishnan I suggest using numerical solution. There is no guarantee that Maple will be able to solve these 9 odes symbolically. I tried, but didn't have the patience to wait long enough to see whether Maple finally gives up.
Numerical solution gives no problem.

@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.