Inserting the matrix on a 1D-input line works. 1D input is also called Maple Notation.
In my view 1D input is way superior to 2D.
This was another example.

@Nia Dutta How did you produce this "approximation"? It is not good anywhere except at one point, basically.

@rlopez Yes, I agree. But if sin(x)/x is the exact solution then that approximation is not only bad, but terrible:
 

numapprox:-infnorm(x^2+423*x^3/3+6789*x^4/4-sin(x)/x,x=0..2);

Value: 28287.54535

@Haqphy12015 Since your a2F(w) values are negative in the beginning the integral will necessarily be negative in the beginning. What to do about that must depend on the physical meaning of a2F(w) (assuming there is one).
Obviously using abs(a2F(w)) instead of a2F(w) will make the integral nonnegative and increasing, but so would e.g. a2F(w)^2.
Which physical meaning do w and a2F(w) have if any?

@tomleslie You are right. A person asking the same question twice (although somewhat differently phrased) at least ought to give us a link to the earlier question and tell us why he wasn't satisfied with the first answer.

@Carl Love Thanks. ApproximateInt fortunately doesn't shy away from symbolic input so it can be used to show your point:
 

restart;
S1:=Student:-Calculus1:-ApproximateInt(f(x),x=a..b,method=right,partition=1);
S2:=Student:-Calculus1:-ApproximateInt(f(x),x=a..b,method=left,partition=1);
factor((S1+S2)/2);
Student:-Calculus1:-ApproximateInt(f(x),x=a..b,method=trapezoid,partition=1);

@acer Thanks for reminding me of the package Interpolation, which is new to Maple 2018.
Also for pointing out that the procedural version of numerical integration is much faster than the nonprocedural that I was using.

@Gabriel samaila Your code resembles code I recall seen written by Carl Love in MaplePrimes, but seems garbled.
If I'm right can you give the link?

@torabi If you go to the help page ?pdsolve/numeric you will see that 2 variables are allowed only, and there is the further restriction that your system consists of time-based partial differential equations.

@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) ?