@bfathi Could you tell me what the problem was?
The code I gave in full should work also in Maple 14. I did it in Maple 2015.

@bfathi OK, but you really don't need this procedure. odeplot will do the job as I described above.

@maple fan As Kitonum pointed out, obviously taking only the first 3 and the last 3 terms of the truncated series for cosine (called p by you) is not very likely to look like cos(70) at all. Why should it?

Since the signs of the cosine series alternate in signs there is no reason to expect s evaluated (one way or the other) at 70 is going to be between -1 and 1.

You mention this as a problem in a book on Mathematica. But you don't mention what the point of the problem is.

Could it be that you are misunderstanding the Mathematica code and that your translation into Maple therefore is incorrect?

Just a question: If you are going to have 10^4 - 10^5 odes, then you are obviously generating them in some programmatic way. Are you varying some parameters or what is the process?

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