@torabi I let the following run for over 10 minutes:
 

pdsolve({seq(PDE[i],i=1..5)});

Then I interrupted it.
Numerical solution in Maple can be done for pdes in 2 variables only.

@Lali_miani In the palette i, j, and I all mean the imaginary unit. If you click on 'i' and your input is 1-D input (aka Maple Input) the result is I, provided that your setting is the default: interface(imaginaryunit=I).
In 2-D input what you get by clicking 'i' is in fact 'i' in the input line, but try to execute the line. The answer will be I, again provided interface(imaginaryunit=I).

MaplePrimes doesn't show this image so I can see it.
Please provide code in any case instead of images.

@acer Sorry, I was too hasty. I just saw the 3.6 minutes.
I was very much aware of the important difference in the order of the ranges pointed out by you in your post in December and was busy finding that post so I could use it.
The other point about not setting the parameter if it hasn't change is not so obvious, but must also be kept in (my) mind.

@acer Did you really mean 3.6 minutes?
On my several years old machine your code takes 10 seconds.
I was just about to answer using your post:
https://www.mapleprimes.com/posts/209925-3D-Plot-Of-Numeric-ODE-With-Parameter
### Since it is much faster to use a real system, I shall bring that version using your notation from that post:
 

restart; # real version
N:=1:M:=sqrt(N*(N+1)):N1:=1+N:w:=10: # f gone
                             
dsys:={diff(z(t),t)=-(N1+M*cos(2*w*t))*z(t)-1+f*(x(t)+y(t)), diff(x(t),t)=-(N1-I*w-2*M*exp(-2*I*w*t))*x(t)-f*(N1+(z(t)))-2*f*M*exp(2*I*w*t),diff(y(t),t)=-(N1+I*w-2*M*exp(2*I*w*t))*y(t)-f*(N1+(z(t)))-2*f*M*exp(-2*I*w*t)};

sys:=subs(x(t)=ux(t)+I*vx(t),y(t)=uy(t)+I*vy(t),z(t)=uz(t)+I*vz(t),dsys);
sysR:=evalc(Re~(sys));
sysI:=evalc(Im~(sys));
SYS:=sysR union sysI;
INI:= {ux(0)=0.5,uy(0)=0.5,uz(0)=0,vx(0)=0,vy(0)=0,vz(0)=0};
res:=dsolve(SYS union INI,numeric,output=listprocedure,parameters=[f]);
Z := eval(uz(t), res);
######################
# acer: https://www.mapleprimes.com/posts/209925-3D-Plot-Of-Numeric-ODE-With-Parameter  (2018)
Q := module()
       local ModuleApply,paramloc;
       ModuleApply := proc(par,var)
         if not (par::numeric and var::numeric) then
           return 'procname'(args);
         end if;
         if paramloc <> par then
           paramloc := par;
           Z(parameters=[f=paramloc]);
         end if;
         Z(var)
       end proc:
end module:
CodeTools:-Usage(plot3d([t,f,Q(f,t)], f=1..10, t=0..3, grid=[29,49],labels=["t","f","z(t)"]));

This takes less than 0.5 seconds on my machine.
I'm still using the old parameters.
You used:

N:=0:M:=0:N1:=1+N:w:=10:

@adel-00 Could you elaborate on what you mean by a 3d plot of g(z(t),t,f) ?
What is g?
And do you mean to vary the parameter f (which had the value 1) ?

@bliengme It appears that the book you are reading has the wrong result (or that you copied it incorrectly).
The notation Q[inf] seems to suggest that that should be the limit of Q(t) as t-> infinity.
That interpretation agrees with the stated result if r > 0.
But the same Q[inf] appears in P: That could be an error in copying from the book (or a misprint in the book).
Consider this where I test if Q__inf is just a constant and is the same appearing in both places:
 

restart;
P := t->Q__inf*r/(2 + 2*cosh(-r*t + b)); # Q[inf]
res:=int(P(s), s = 0 .. t); 
res1:=convert(res,exp);
res0:=Q__inf/(1+exp(b-r*t));
# A simple test:
eval(res1-res0,{r=7,b=3,t=5,Q__inf=1});
evalf[30](%);
# Numerical integration:
res2:=Int(P(s), s = 0 .. t); 
evalf[30](eval(res2-res0,{r=7,b=3,t=5,Q__inf=1}));

The results show that the book (or your copying) is wrong.
## NOTE: The two res0 and res differ by a constant since their derivatives are equal:
 

res0:=Q__inf/(1+exp(b-r*t));
##
diff(res0,t)-P(t); 
convert(%,exp);  # 0

@torabi SYSPDE is a set having 20 members. So far so good.
But take a look at e.g. SYSPDE[20].
I get:
-9.585649527*10^13*xi[3, 3] + 1.171579387*10^14*xi[3, 4]
     + 1.437847430*10^14*xi[4, 3] - 1.757369082*10^14*xi[4, 4]
     - 2.955249835*10^16*kappa[3, 3] + 3.149063384*10^15*kappa[3, 4]
     + 1.235128694*10^16*kappa[4, 3] - 4.260288677*10^13*kappa[4, 4] + Matrix(
    6, 6, [[-0.5521334135*10^15, 0., 0., 0., 0., 0.], [0., -0.5521334135*10^15
     - 2532.446978*kappa[3, 3] + 1266.223489*kappa[3, 4]
     + 1266.223489*kappa[4, 3] - 633.1117444*kappa[4, 4],
    -2849.002850*kappa[3, 3] + 1424.501425*kappa[4, 3],
    -2849.002850*kappa[3, 4] + 1424.501425*kappa[4, 4], 1266.223489*kappa[3, 3]
     - 2532.446978*kappa[3, 4] - 633.1117444*kappa[4, 3]
     + 1266.223489*kappa[4, 4], 0.], [0.,
    -2849.002850*kappa[3, 3] + 1424.501425*kappa[3, 4],
    -0.5521334135*10^15 - 3205.128206*kappa[3, 3], -3205.128206*kappa[3, 4],
    1424.501425*kappa[3, 3] - 2849.002850*kappa[3, 4], 0.], [0.,
    -2849.002850*kappa[4, 3] + 1424.501425*kappa[4, 4],
    -3205.128206*kappa[4, 3], -0.5521334135*10^15 - 3205.128206*kappa[4, 4],
    1424.501425*kappa[4, 3] - 2849.002850*kappa[4, 4], 0.], [0.,
    1266.223489*kappa[3, 3] - 633.1117444*kappa[3, 4]
     - 2532.446978*kappa[4, 3] + 1266.223489*kappa[4, 4],
    1424.501425*kappa[3, 3] - 2849.002850*kappa[4, 3],
    1424.501425*kappa[3, 4] - 2849.002850*kappa[4, 4], -0.5521334135*10^15
     - 633.1117444*kappa[3, 3] + 1266.223489*kappa[3, 4]
     + 1266.223489*kappa[4, 3] - 2532.446978*kappa[4, 4], 0.],
    [0., 0., 0., 0., 0., -0.5521334135*10^15]]) + 1.022469284*10^15*w[4, 3]
     - 1.533703927*10^15*w[4, 4] - 5.410566626*10^15*w[3, 3]
     + 7.157284986*10^15*w[3, 4]

You will notice that this is a sum of algebraic terms and a Matrix.
Try:
 

for i to nops(SYSPDE[20]) do i = op(i,SYSPDE[20]) end do;

whattype(op(9,SYSPDE[20]));

There may be similar problems for other SYSPDE[i].

@Dark Energy To plot diff(a(t),t)/a(t) use that diff(a(t),t) is given as the right hand side of ode1:
 

## Continuing with ode1 as defined previously:
ode1;
res:=dsolve({ode1,a(0)=1},numeric,parameters=[c,w,z]);
res(parameters=[1,2,3]); # Setting parameters of choice
plots:-odeplot(res,[t,a(t)^2-1],-1..10); # The easy one
## For the next you have to evaluate the right hand side at the parameter values:
plots:-odeplot(res,[t,eval(rhs(ode1),{c=1,w=2,z=3})/a(t)],0..10,thickness=3,labels=[t,diff(a(t),t)/a(t)]);

Note about the somewhat broken line from ca. 6 to 10:
rhs(ode1) is zero here, but this is a numerical computation, so roundoff happens. Thus some of the t values could result in imaginary numbers because of the square root. They can't be plotted, thus the holes.
 

eval(rhs(ode1),{c=1,w=2,z=3});
eval(%,res(8)); # An example where I get : 2.76925989236618*10^(-8)*I

In the help page for pdsolve under Parameters we find about 'build':
"(optional) try to build an explicit expression for the indeterminate function, no matter what the generality of the solution found"
(my emphasis).

After just a glance at table 3 (thanks to Tom Leslie) I find that what it says is:
"Proposed approach (approximate, Maple)".
Another Maple entry is:
"Proposed approach (Maple)".
Obviously, by "proposed approach" is meant methods proposed by the authors. So the difference between the two mentioned must be found by reading the article, which I haven't done.
There is no approximate Maple as opposed to just plain Maple.

@tomleslie

I agree with your three points:

1. Put the restart statement in its own execution group
2. Don't use 2D input
3. Don't use document mode

@Carl Love Even in document mode and with 2d input I get -1 as I should.

@Jodelkoenig I just tried in Maple 2017.3 with restart separated as Carl describes. I had no problems.
Which Maple 2017 do you have?
Mine is:
interface(version);
`Standard Worksheet Interface, Maple 2017.3, Windows 10, September 27 2017 Build ID 1265877`
and

kernelopts(version);
`Maple 2017.3, X86 64 WINDOWS, Sep 27 2017, Build ID 1265877`

@awass There is no file attached.