@In-Jee Jeong Since it is not clear what you are actually doing I think you ought to present us with the details, perhaps as an uploaded worksheet.

It is not clear what you are actually doing because in the question itself the word 'numeric' does not appear in the dsolve command as tomleslie points out, and yet you are using odeplot there.
Secondly, in your reply to bogo you write plot(b[0](t),t=0..1). For a numerical solution that will only work if you are using output=listprocedure and b[0] is assigned one of the solution procedures.

@tomleslie I agree. The problem is (again) the 2D input, where spacing is interpreted as multiplication, which apparently is indeed intended in this case. But then there are (it appears) two cases of missing a space after PB.

Since everybody knows that if an expession F is an antiderivative of the expression f on some interval I then so is F+C, where C is any constant, there is no need for Maple to add a constant. What would be the purpose? Knowing one antiderivative (on some interval) means knowing all.
You just have to know that Maple doesn't add such a constant.

I noticed that in the help page for int in Maple 2015 this is pointed out explicitly.

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