@trottaleonardo OK, so you tried a different approach. But your code makes a 1 x nops(Z) matrix.
Revise it to (assuming w is the same as your Z):

convert(w,list);
`[]`~(%); #To turn the list into a listlist
convert(%,Matrix);

In the example you gave w was a sum of terms (type `+`). I'm assuming this is also the case with Z. If nops(Z) = 2 and Z is of type `+` then there are only 2 terms. So you get only two lines.

By the way, did you actually try my first code on your own example?

Your code exactly as written (except for supplying a missing semicolon at the end) worked for me.  I'm using ghostscript and ghostview for viewing. 

Try

Digits:=25;

and a much smaller range (0..0.0001).

Added later:

In order to see that X is not constant (although it appears to be at first sight) you can plot like this:

rg:=0..0.0001;
X0:=10^20:
init:=R(0)=0,X(0)=X0:
res:=dsolve({sys,init},[X(x),R(x)],numeric,range=rg);
plots:-odeplot(res,[x,R(x)],rg);
plots:-odeplot(res,[x,ln(X(x)/X0)],rg,thickness=4,view=-1.5e-17..0);

#Alternatively, make a simple change of variableX(x) = X0 + xi(x):

X0:='X0':
sys2:=PDEtools:-dchange({X(x)=X0+xi(x)},{sys},[xi]);
Digits:=25:
X0:=10^20:
init2:=xi(0)=0,R(0)=0:
rg:=0..0.0001;
res:=dsolve(sys2 union {init2},[xi(x),R(x)],numeric,range=rg);
plots:-odeplot(res,[x,R(x)],rg);
plots:-odeplot(res,[x,xi(x)],rg,thickness=4);

Try

Digits:=25;

and a much smaller range (0..0.0001).

Added later:

In order to see that X is not constant (although it appears to be at first sight) you can plot like this:

rg:=0..0.0001;
X0:=10^20:
init:=R(0)=0,X(0)=X0:
res:=dsolve({sys,init},[X(x),R(x)],numeric,range=rg);
plots:-odeplot(res,[x,R(x)],rg);
plots:-odeplot(res,[x,ln(X(x)/X0)],rg,thickness=4,view=-1.5e-17..0);

#Alternatively, make a simple change of variableX(x) = X0 + xi(x):

X0:='X0':
sys2:=PDEtools:-dchange({X(x)=X0+xi(x)},{sys},[xi]);
Digits:=25:
X0:=10^20:
init2:=xi(0)=0,R(0)=0:
rg:=0..0.0001;
res:=dsolve(sys2 union {init2},[xi(x),R(x)],numeric,range=rg);
plots:-odeplot(res,[x,R(x)],rg);
plots:-odeplot(res,[x,xi(x)],rg,thickness=4);

@Markiyan Hirnyk Certainly bug fixing should be done by employees of Maple. But while we are waiting for Maple 16.03 we will have to live with what we got.
I'm afraid that it is unlikely that any update will appear before Maple 17 is released (and probably not after either).
If that is so, people who don't buy Maple 17 will have to live with that bug forever.

Maybe you should have waited till you could upload the mw-file.

In one of the links given by acer:

http://www.mapleprimes.com/questions/142322-Why-Has-The-Plot-Command-Deteriorated#comment142326

I suggest a more general workaround:

R:=proc(p) local f;
   f:=proc(curve)
      evalindets(curve,Or(listlist,Matrix),m->ListTools:-Split(hastype,convert(m,listlist),undefined))
   end proc;
   evalindets(p,specfunc(anything,CURVES),f);
end proc;

R(A); #A being the plot

Did you try copy and paste of the lines I wrote? Why does the error message mention an F? The unknown function has the name X.

Did you try copy and paste of the lines I wrote? Why does the error message mention an F? The unknown function has the name X.

@Simon Willerton Using indexed names sometimes gives (unexpected) problems because of the meaning an attached [..] also has.

Try

LinearMultivariateSystem({ a[1] + a[2] = 1, a[1] > 0, a[2] > 1, a[3] > 0 }, {a[1],a[2],a[3]});

and you will see a weird error message.

To construct names like f1, f2, etc. you can use cat, as in

seq(cat(f,k),k=1..3);

Or after the fact:

evalindets([{ a[1] + a[2] = 1, a[1] > 0, a[2] > 1, a[3] > 0 }, {a[1],a[2],a[3]}],indexed,x->cat(op(0,x),op(1,x)));

and then

LinearMultivariateSystem(op(%));

@Simon Willerton Using indexed names sometimes gives (unexpected) problems because of the meaning an attached [..] also has.

Try

LinearMultivariateSystem({ a[1] + a[2] = 1, a[1] > 0, a[2] > 1, a[3] > 0 }, {a[1],a[2],a[3]});

and you will see a weird error message.

To construct names like f1, f2, etc. you can use cat, as in

seq(cat(f,k),k=1..3);

Or after the fact:

evalindets([{ a[1] + a[2] = 1, a[1] > 0, a[2] > 1, a[3] > 0 }, {a[1],a[2],a[3]}],indexed,x->cat(op(0,x),op(1,x)));

and then

LinearMultivariateSystem(op(%));

In Maple 15 and 16 the result is the correct NULL (i.e. no solution):

solve({ a + b = 1, a > 0, b > 1 }, {a,b});

However, the result of

solve({ a + b + c = 1, a > 0, b > 1, c > 0 }, {a,b,c});

is also wrong.

It is the same in Maple 16.

As always you are more likely to get useful answers if you provide us with all the details.

Could you give us the pde and the ibc the way you gave it to Maple? (Text please).