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

Unintentionally, I had a range in both of the loops above. It has now been removed in the first loop.
Also notice that the result of maximize(cos(t*a),t=1..20) is wrong when a is just a name.

@Markiyan Hirnyk Yes, that must be a bug:

for k in [Pi,gamma,sqrt(2),sin(1),123,0.4567,a,b] do
  try
   res:=maximize(cos(t*k)); #Revised: Here was a range, now no range
   print(res);
  catch:
   print(k,lastexception[2]);
  end try;
end do;
Error for Pi,gamma,sqrt(2),sin(1).
##Better results with a range:
for k in [Pi,gamma,sqrt(2),sin(1),123,0.4567,a,b] do
  try
   res:=maximize(cos(t*k),t=1..20);
   print(res);
  catch:
   print(k,lastexception[2]);
  end try;
end do;
## The a (and b) results are wrong, though:
maximize(cos(t*a),t=1..20);#Result:  max(cos(a), cos(20 a))
eval(%,a=Pi/9); #Not equal to 1 although
maximize(cos(t*Pi/9),t=1..20); #Result 1, which is correct.

@vv Looks like a bug to me. I shall submit an SCR.

It helps to read the help page for DEtools[DEplot].

But first of all: What is a(t)?
Why do you use x and y when the dependent variables are called b and c?
Where is the list of initial conditions?

Try (non-numerically):

dsolve(sys,{b(t),c(t)});

@vv You say it used to work that a and a(t) are taken to be the same after the alias.
That would be quite weird. Consider a(t)(t)(t)(t). Should that then reduce to a?
In fact it doesn't work like that in Maple 15. I just checked.

In Maple 15 and Maple 2015.2:
restart;
alias(a=a(t));
whattype(a); #function
whattype(op(0,a(t))); #function
eval(a(t),t=5); #a(5)(5)

@Bendesarts For the fun of it I tried the following, where the fun part is in evalindets.
You can use any unassigned name instead of %f, e.g. just f.

restart;
eq1 := l1*cos(theta(t))+l2*sin(beta(t))-x(t);
eq2 := l1*sin(theta(t))-l2*cos(beta(t));

T,S:=op(map[3](evalindets,[eq1,eq2],function(identical(t)),[s->op(0,s),s->%f(op(0,s)=s)]));
S:=indets(S,`=`);
Phi:=subs(S,VectorCalculus[Jacobian](T,[theta,x,beta]));