@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]));

@max125 The manuals are good. You can find them here:

http://www.maplesoft.com/documentation_center/

@max125 When you do (in Maple input):

g(x):=x^2;

you are (implicitly) creating a procedure g with option remember and assigning to its memory table:

eval(g);
            
Notice that the body of this procedure i.e. just 'procname(args)' makes this procedure return unevaluated if the input is not in the memory table of g. The assignment you made has been saved in the memory table:

op(4,eval(g));
               

Here x is the index and x^2 is the entry.
Try one more assignment:
g(87):=56;
op(4,eval(g));
                  
To define a function g given by g(x)=x^2 for all x to:
g:=x->x^2;
# or use the more explicit version:
g:=proc(x) x^2 end proc;

@adel-00 The amplitude of the oscillations is very small. I commented on that in my answer above, where I ended the Maple code with the comment ##No perceptible wiggles.

Try replacing 0.9..1 by 0.99999..1 (as I did in my comment to vv above). Then you get:

@vv Yes using
dsol1(parameters=[10^6]);
plots:-odeplot(dsol1,[t,u(t)],0.99999..1,numpoints=10000);
## I get:


which seems to have the same period as sin(2*10^6*t), i.e. Pi*10^(-6) :
plot(sin(2*10^6*t),t=0.99999..1);

@adel-00 Answers to your questions:

1) the size option only works in recent versions of Maple. Just remove it.

2) Nothing wrong with setting omega=10^6 at the top, but then you loose the opportunity to analyze numerically what the effect of that extremely rapid oscillation has.
If you don't want to experiment with omega you should just use:
dsol1 :=dsolve({dsys1,ini1},numeric, abserr=1e-9, relerr=1e-8,maxfun=0);
plots:-odeplot(dsol1,[[t,u(t)],[t,v(t)],[t,w(t)]],0..1,numpoints=10000);

3) 20 minutes doesn't sound wrong. On my computer with Maple 2015.2 the last command took 48 seconds. You have an earlier Maple version and you may have a slower computer.