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.

@Carl Love Yes, I wouldn't have thought of that, but isolating
 zip(solve, [(x-2)/3 , (y-1)/4 , (z-3)/3] =~ t, [x,y,z]);
made me understand.

@vv Certainly, yes. If this it not desirable then the first method should be chosen, i.e.

f:=proc(theta) proc(x) theta*x end proc end proc;

You have an initial value problem for the 1D wave equation. The solution is given by the result of
pdsolve({Eq1,u(x,0)=f(x),D[2](u)(x,0)=g(x)});

This is THE solution. So what more do you want?

@maple I edited my original answer above to include a shooting method. Take a look at that.
The conclusion I draw from that is that the original bvp solution is pretty good.