@Carl Love method=lsode[adamsfull] seems to work fine:

SollsodeFull := dsolve({ODEs}, numeric, method=lsode[adamsfull]);

SollsodeFull(0.5);

returns

[x = .5, f[0, 0](x) = 0.453528451308833e-2, f[0, 1](x) = -11.5167207961770]

@patient An immediate observation about your worksheet res_syst.
There are two occurrences of v where these should be v(z).

Thus the odes should be:

odes:=0.4e-2*z*diff(u(z),z,z)+(z-8*z*diff(v(z), z)+0.6e-2)*diff(u(z), z)+(3/2-12*diff(v(z), z)-8*z*diff(v(z),z,z))*u(z) = 0, .8*z*(diff(v(z), z)^2+v(z)*diff(v(z),z,z))+(2*(.6+z))*v(z)*diff(v(z), z)+v(z)*u(z)+v(z)^2 = 0;

With that system and initial conditions
ics1:=D(u)(eps) = 0, D(v)(eps) = 0, u(eps) = 1, v(eps) = 1;

with eps=0 the system has a singularity at 0.

You can play with eps=0.001 or similar.

@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=%);