@Kitonum For the integral int(abs(cos(n*x)), x=0..Pi) try
restart;
int(abs(cos(n*x)), x=0..Pi,method=_RETURNVERBOSE);
#All methods fail.
So it must be an `assuming` problem.
I have always been suspicious of assuming equality.
Try
restart;
assume(n=0);
int(abs(cos(n*x)), x=0..Pi,method=_RETURNVERBOSE); # ftoc = 0, ftocms = 0, all other fail.
about(n);
B:=int(abs(cos(n*x)),x); # Result: sin(n*x)/n
is(n=0); #true
eval(B,x=Pi/2); # 0
restart;
assume(n::integer);
int(abs(cos(n*x)), x=0..Pi,method=_RETURNVERBOSE); #All fail
about(n);
int(abs(cos(n*x)),x); #signum(cos(n*x))*sin(n*x)/n
################
## Totally apart from integrals, take a look at the following, where my question is: What kind of answer do we expect or want?
restart;
assume(n=0);
n^2; # We are getting n^2
eval(%); # no change
sin(n*Pi); # 0
## Why the difference? Is it reasonable?

@awass

Well, method = laplace is not the default.
Notice that when using method=laplace and type=series (or just series) you need to give the unknown function(s) as the second argument:
ode:=diff(x(t),t)=x(t);
dsolve(ode,x(t),method=laplace); #Notice x(0) used as the arbitrary constant (natural with laplace)
dsolve(ode,x(t),series,order=4); #Notice x(0) used as the arbitrary constant (since t=0 is used for expansion point)
dsolve(ode); #No need for x(t), but OK with it.
#_C1 used as arbitrary constant as it also is when the ode is nonlinear (consistency!)
dsolve({ode,x(0)=1},numeric); #No need for x(t), but OK with it
I understand from Edgardo Cheb-Terrab that the difference wrt. the need for x(t) is simply due to the fact that different people are in charge of these things.

@Markiyan Hirnyk I wanted to get rid of the square root, so needed to solve sqrt(x+ln(x))=u for x.

@Rouben Rostamian  In Maple 2016 the code executes without error. In Maple 2015.2 we need the assumption eps>0 as vv mentions. But see note below.
An interesting change from Maple 2015.2 to Maple 2016.
You could say that as long as the integrals are inert there shouldn't be any complaining, the burden of validity rests on the user, as it does for several operations in Maple.
Or you could say that if the one integral is convergent (exists) then so is (does) the other.
##
Note. The burden rests on the user:
In Maple 2015.2 this works too:
Change(J1,s=-t) assuming eps<0;
In this case the integral J1 is divergent, but so is the integral emerging from Change.

@Declan Matrices are indexed from 1 and up. Arrays are more flexible. So you could do:

f := (i,j) -> `if`(j=1, i^2, i^3):
A:=Array(5..14,1..2,f);
A[5,..]; # .. means all of the second range i.e. in this case 1..2
A[7,2];

@adel-00 Your expression w1*(-1/3) is imaginary, as you also seem to realize since you are taken its conjugate in defining f:
f:=conjugate(w1)*(-1/3):

contourplot deals with real-valued expressions only. Thus you can use it on Re(f) or on Im(f) or on abs(f) as you please.
You are asking for only one contour level, -1. If you use Re(f) that level is not reached in the rectangle you chose.

@adel-00 The solutions to the ode for z are oscillatory.

We can prove that there exists exactly one solution which is periodic with period Pi/omega (i.e. the period of cos(2*omega*t) ).
I believe that we should be able to prove that any other solution will approach that solution as t->infinity.
((Easy proof added at the bottom))
Thus limit( z(t), t=infinity) doesn't exist, no matter which solution z(t) you have.
The amplitude of the oscillations is, however, very small when omega=10^6 as you have, and the oscillations are indeed about a median of -1/3 (rather closely).

restart;
ode:=diff(z(t),t)=-(N1+M*cos(2*omega*t))*z(t)-1;
res:=dsolve({ode,z(0)=z0});
cond:=eval(res,{z(t)=z0,t=Pi/omega}); #Condition on z0 for periodicity with period Pi/omega
zp:=solve(cond,z0); #Exactly one solution for z0.
##Rewriting the integral:
zp:=IntegrationTools:-Change(zp,_z1=t/2/omega,[t]);
expand(zp);
normal(%);
zp:=combine(%,exp);
#Now with your parameters:
zp0:=evalf[35](eval(zp,{omega=10^6,N1=3,M=sqrt(2)}));
evalf[25](eval(zp,{omega=10^6,N1=3,M=sqrt(2)}));
Digits:=25:
resnum:=dsolve({eval(ode,{omega=10^6,N1=3,M=sqrt(2)}),z(0)=zp0},numeric,abserr=1e-12,relerr=1e-12,output=listprocedure);
plots:-odeplot(resnum,[t,z(t)],0..Pi*10^(-6));
Z:=subs(resnum,z(t)):
mx:=Optimization:-Maximize(Z,0..Pi*10^(-6));
mn:=Optimization:-Minimize(Z,0..Pi*10^(-6));
mx[1]-mn[1]; #Twice the amplitude of the oscillations: 5.94*10^(-7)
(mx[1]+mn[1])/2; # Oscillates about that value: -.3333332720294223127799336
##So -1/3 seems close enough.
################################
##Proof that any solution approaches the periodic solution determined above:
The difference between any two solutions satisfies the homogeneous equation, which is easily solved:
odeH:=diff(z(t), t) = -(N1+M*cos(2*omega*t))*z(t);
resH:=dsolve({odeH,z(0)=z0});
limit(rhs(resH),t=infinity) assuming omega>0,N1>0; #Result 0
#We have an upper bound for abs(z(t)):
eval(rhs(resH),{sin=-1,z0=abs(z0)});
combine(expand(%));
## result   abs(z0)*exp(-N1*t+(1/2)*M/omega)

@adel-00 Your system dsys is not an autonomous system. Its right hand sides depend on t explicitly, i.e. not only implicitly through x(t), y(t), and z(t), but explicitly through exp(2*I*omega*t) and cos(2*omega*t).
So res:=solve(eqs,{x,y,z}); doesn't have any significance that I can see.

Your ode for z doesn't depend on x and y so can be solved independently of those by dsolve
The problem is that the integral in the result cannot be found:
restart;
ode:=diff(z(t),t)=-(N1+M*cos(2*omega*t))*z(t)-1;
dsolve({ode,z(0)=z0});
We see that we may set z0=0 if we are only interested in the limit of z(t) as t->infinity

@Jenser As far as the derivative sort goes, the definition of S has to reflect the number of derivatives in L, so the number 4 has to be replaced by the more general nops(L) or numelems(L):
S:=seq(res[ListTools:-Search(L[i],[k])]*L[i],i=1..nops(L)),res[ListTools:-Search(1,[k])];
Then the code should work.
##
As for indices as in a__1 you can do:
ex:= a__2+a__3+a__1;
sort(ex,order=tdeg(a__3,a__2,a__1)); #Assuming you want this done descendingly
#or by constructing the sequence used in tdeg() first:
S:=seq(cat(a__,i),i=3..1,-1);
sort(ex,order=tdeg(S));
##
As for improvements to the document mode in Maple 2016 I cannot tell you anything since I never use Document mode.
##
Finally to your comment about sorting not being unusual.
Mathematically, there is no difference between a+b and b+a. Maple keeps a copy of only one of these in memory for efficiency reasons. As long as we are talking mathematics there is no reason to waste time on sorting.
For presentation use sorting makes sense, though.

@Kitonum I think this is the time to point out that devices like the ones we are talking about may be OK at the interactive level, where you can see immediately what you are getting. In a programmatic context such devices should not be relied upon.
A couple of days ago I reviewed a worksheet in which (at the interactive level) I referred to something like op([1,1,1],A).
Because A in Maple 2016 looked slightly different op([1,1,1],A) no longer referred to what I wanted, but I could see that.

The ordering of sets in older versions was quite unpredictable. A change was made (maybe starting in Maple 17?) to make the ordering of sets consistent from session to session.
Note.
On a 32 bit computer in the standard interface I tried the scheme on Maple 15, 16, 17, 18.02, and 2015.2.
In Maple 15 and 16 I get from
L:=convert(indets(Expr,specfunc(anything,diff)),list);

 whereas in Maple 17, 18, and 2015 I get

(in addition in Maple 17 and later you don't need the 'anything' in specfunc).
The scheme works in 15 and 16 if you make use of the different ordering of L:
S:=seq(res[ListTools:-Search(L[i],[k])]*L[i],i=4..1,-1),res[ListTools:-Search(1,[k])];
In Maple 17 and 2015 the scheme works unaltered, but it doesn't work in  Maple 18.02 (well sometimes).

@Kitonum My solution does work in Standard Worksheet Maple 2015.2 as well as in Maple 2016 (at least on my computer). I just double checked.

On this computer I don't have access to the classic interface.

But I agree, your latest solution doesn't work in Maple 2016 Standard interface nor in Maple 2015.2 Standard interface.

The problem is that S1 is not in the order we want.

By using indets(Expr,specfunc(diff)) I get the lexical ordering: diff(diff(x(t),t),t) will precede diff(diff(y(t),t),t), which will precede diff(x(t),t), which again will precede diff(y(t),t).

In your case the ordering of S1 is:
[a*(diff(x(t), t)), r*(diff(x(t), t, t)), a*(diff(y(t), t)), b*(diff(y(t), t, t)), b*x(t), c*y(t)]
which isn't lexical (for whatever reason). If you replace S1 by the set version
S1:=op~({selectremove(s->has(s, diff), [op(Expr)])});

then the ordering is lexical:
{a*(diff(x(t), t)), a*(diff(y(t), t)), b*(diff(y(t), t, t)), b*x(t), c*y(t), r*(diff(x(t), t, t))}
but not what you want. The problem is the coefficients a, b, c, r.

P.S. Let me add that I always use the Standard interface in worksheet mode and Maple input (aka 1D math input).

@Kitonum I tried a different version:

restart;
Expr:=a*diff(x(t),t)+b*x(t)+r*diff(x(t),t,t)+a*diff(y(t),t)+b*diff(y(t),t,t)+c*y(t);
L:=convert(indets(Expr,specfunc(diff)),list);
res:=coeffs(Expr,L,k);
k;
S:=seq(res[ListTools:-Search(L[i],[k])]*L[i],i=1..4),res[ListTools:-Search(1,[k])];
## Now observe the difference between these 3 versions of the ode with `` in different places:
add(``(i),i=S); #OK, but parentheses around all terms of course.
``(add(i,i=S[1..4]))+S[-1]; #Remarkable that the order of S is kept inside the `` function as in your version.
``(add(i,i=S)); # Order lost inside the `` function.

Comment. The reason for the loss of the order in the last version (and in the 4th version add(i,i=S) ) must be simply that the original expression Expr is recognized as being already in memory and so that is used.
The reason then that sorting the indexed version works must be that  a[2]+a[3]+a[1] is a polynomial in 3 variables for which sorting is possible (and is a destructive operation, i.e. is done in place).

@Jenser You can revise Carl's lines like this:

restart;
ex:= a[2]+a[3]+a[1];
sort(ex, sort([op(ex)], key=-op));

@Declan Probably there is a way since the caption in this plot prints nicely:

plot(sin,0..Pi,caption=typeset("The solution to ", a^x=b, " is ",x=1323));

I hope somebody more knowledgeable in these matters than me will read this and help you.

@tomleslie I'm assuming that your comment was meant for vv. In my understanding of his question

"Maple 2015 can compute it, but with a little help.
What about Maple 2016?"

both Markiyan and I answered that Maple 2016 can do the integral with a little help as can Maple 2015.