@vv I'm not sure what you mean by equality not being the real issue.
Certainly as you say we have from
assume(p<=0,p>=0);
about(p);
the response
Originally p, renamed p~:
  is assumed to be: 0

So we are actually assuming that p = 0, right?
I'm not saying that it is wrong for assuming to accept equality, I'm just asking what the designers want the system to respond to the fact that the assumptions imply equality.
As it is now the response is not easy to guess.

Should anybody return to this let me give the last revision that I have

restart;
assumingoverload:=subsop(3=overload,eval(`assuming`)):
st:=convert(eval(assumingoverload),string):
st2:=StringTools:-Substitute( st, st[1..19], "proc(x::uneval,a::list(Not(satisfies(x->hastype(x,name=constant)))))" ):
STR2:=" if type(a,'list') then assumptions := op(a)":
st3:=StringTools:-Substitute(st2,STR2,cat(STR2[1..-3],"ListTools:-Flatten(evalindets(a,set,[op])))")):
assumingoverload:=parse(st3):

assumingoverload:=parse(st3):
`assuming`:=overload(
       [
          proc(x::uneval,par1::And(list,satisfies(x->hastype(x,name=constant)))) option overload;
           local par,par2,res;
           par:=ListTools:-Flatten(evalindets(par1,set,[op]));
           par,par2:=selectremove(hastype,par,name=constant);
           assign(par);
           try
            res:=eval(x) assuming op(par2);
           catch:
            error;
           finally  
            map(unassign@lhs,eval(par,2));
           end try;
           op(eval(res,par)) assuming op(par2);
          end proc,
          eval(assumingoverload)
        ]
):

#################
Examples.
plot(a*sin(x),x=0..b) assuming {a=-2,b=7};
param := {a=2,b=7}:
plot(a*sin(x),x=0..b) assuming param;
plot(a*sin(x),x=0..b) assuming a=-2,b=7;
plot(a*sin(x),x=0..b) assuming [a=-2],{b=Pi};
simplify(sqrt((x*y)^2)) assuming x>0,y>0;
simplify(sqrt((x*y)^2)) assuming x=10,y>0;
int(exp(a*t),t=0..b) assuming a<0,b=infinity;
int(exp(b*a*t),t=0..infinity) assuming b=2,a<0;
int(exp(b*a*t),t=0..infinity) assuming {a<0,b>2};
A:=Matrix([[a,b],[c,d]]);
LinearAlgebra:-Eigenvectors(A) assuming d=0;

@Markiyan Hirnyk Thank you Markiyan. Yes that's the one.

@vv Several years ago I wrote contrived an overloading of assuming to work as eval when assuming equality.

That works fine, but the question is what the designers of `assuming` and assume intend assuming equality to mean.

I will formulate a question in MaplePrimes with the title "What does assume n=7 mean?"

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