You could try following the advice given, i.e. using optional arguments
maxmesh=256 (or higher, 8192 being the highest available) and/or abserr=1e-3 (just an example).
That may or may not help. If not, let us see the whole problem.

@Vesnog Here is a simple example:

eq:=tan(x)=a*x;
sol:=[seq(fsolve(eq,x=0..Pi/2,avoid={x=0}),a=[1.1,1.6,Pi])];

Your expression RHS-LHS is not a polynomial in x, thus fsolve will not return all 3 solutions on Pi/2..lim.
The method of splitting you are using works fine, but obiously has the disadvantage of relying on the graph or at least on knowledge of the singularities of RHS-LHS.

@edrobina I have no problem with your worksheet in Maple 18.

So apparently the set  {a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p} does not only have (unassigned) names a,b,c,...
What is actually in {a, b,  c, d, e, f, g, h, i, j, k, l, m, n, o, p} ?

I think this was a display prpblem in Maple 16 (later corrected as Markiyan points out).
To attempt a verification of this you may try:
taylor(1/(1+z^2), z = 1, 4);
convert(%,polynom);

That it is a display problem is indicated by the fact that the term 1/2 and the term -1 in the incorrect looking result would automatically be simplified to -1/2.

@ewas OK, but you still need to find the four extra boundary conditions relevant to your problem.

There are som missing parameter values. In particular B1 and C1 since they are rather crucial being the coeffcients of highest order. Thus it makes an enormous difference if they are zero or not. On the other hand the information that C2=0 is irrelevant since it doesn't appear in your system.
The spelling is 'dsolve', not 'Dsolve'.
What are the values of az, aw (you wrote 'and'), and ag?
Two of the boundary conditions could be
az * Tz (0) + aw * Tw (0) + ag * Tg (0)=-10, az * Tz (20) + aw * Tw (20) + ag * Tg (20)=10.
You need 4 additional ones.
I don't understand the do-loop.

Since tan(x)  appears in your system and has a singularity at all odd multiples of Pi/2 you should reconsider using the interval x=0..20.

@maple fan I had no problem plotting the solution on the range 0..8.
The plot was rather boring: Very soon the graph became linear with slope 0.1.
What did you expect should happen?
What has mu got to do with the problem?

The following code does the same thing:
restart;
u:=1003-1000*x(t)-30*diff(x(t), t)-25*signum(diff(x(t), t)-.1)-.3*signum(diff(x(t), t))*exp(-2*abs(diff(x(t), t)));
sys := diff(x(t), t, t) = (1-b(t))*u, x(0)=1, D(x)(0) = 0;
stick := [[diff(x(t), t) = .1, abs(1-x(t))< 0.1],b(t)=1];
slip := [[0, abs(1-x(t))>0.1], b(t) = 0];
sol:=dsolve({sys,b(0)=0},numeric,discrete_variables=[b(t)::boolean],events=[stick,slip]);
plots:-odeplot(sol,[t,x(t)],0..0.01);
plots:-odeplot(sol,[t,x(t)],0..8);
plots:-odeplot(sol,[t,diff(x(t),t)-0.1],0.01..8);

@maple fan This behavior is not a bug.
Take F1.
piecewise first examines if t>=0. If that is the case then it returns 0. It does not go on to the next condition to see if that is also satisfied. It only keeps going if previous conditions are not satisfied.
Thus the plot is a plot of the zero function.
If you want the same result as F2 without using 'and' you could do:
F1a:=piecewise(t<2,0,t<2.1, t^3, exp(t));

@maple fan Here is an answer to the first question.
I use two discrete variables, one for turning on and off and one for getting the time where y(t)=30.
restart;
F1:=y(t)^2;
ode:=diff(y(t),t$2)=F1-b1(t)*piecewise(t<b2(t)+0.1,y(t)^3,exp(y(t)));
ic:=y(0)=1.2,D(y)(0)=0;
res:=dsolve({ode,ic,b1(0)=0,b2(0)=0},numeric,discrete_variables=[b1(t)::boolean,b2(t)::float],events=[[y(t)=30,[b1(t)=1,b2(t)=t]]]);
plots:-odeplot(res,[t,y(t)],0..5);
res(3);
plots:-odeplot(res,[t,b2(t)],0..3);
#############
For the second question I took another example for F1 and introduced a counter b3(t).
There are two events. The first is [y(t)=y0,b3(t)=b3(t)+1], which updates b3(t) when y(t)=y0.
The second is conditional and has a double action: [[y(t)=y0,b3(t)=2],[b1(t)=1,b2(t)=t]].
This latter triggers when y(t)=y0 provided b3(t)=2. If so b1(t) is set to 1 and b2(t) to t.

restart;
F1:=sin(y(t)^2);
ode:=diff(y(t),t$2)=F1-b1(t)*piecewise(t<b2(t)+0.1,y(t)^3,exp(y(t)));
ic:=y(0)=1.2,D(y)(0)=0;
y0:=1.9:
res:=dsolve({ode,ic,b1(0)=0,b2(0)=0,b3(0)=0},numeric,discrete_variables=[b1(t)::boolean,b2(t)::float,b3(t)::integer],
  events=[[y(t)=y0,b3(t)=b3(t)+1],[[y(t)=y0,b3(t)=2],[b1(t)=1,b2(t)=t]]]);
plots:-odeplot(res,[[t,y(t)],[t,y0],[t,y0*b1(t)]],0..8,linestyle=[1,2,1]);
res(4);
plots:-odeplot(res,[t,b3(t)],0..4,thickness=3);

@mahmood180 As I wrote you need the output to be
Array(1..N, 0..M , b)
not what you had (which involved b[m,n]).
Secondly, change the assignment
OB(i,M):=...
to
OB[i,M]:=...

@one_man In copying quotes ``were left out of `minus`.

Converting to set and then to list may change the order:
restart;
a := [y+1, x+2, x+3, x+4]; #Notice the y
a := convert(a, set);
a :=`minus`(a, {x+2});
a := convert(a, list); #Notice the order
###########
## Another example:
restart;
a := [x+5, x+2, x+3, x+4];
a := convert(a, set);
a :=a minus {x+2};
a := convert(a, list);

@sirmotion If you follow the link I gave you will see that 'res' refers to the d'Alembert formula.

@ESchiltkamp I have edited your last worksheet to include the numerical approach.
Notice several points:
1. I use the parameters option for dsolve/numeric, see
?dsolve,interactive,parameters
2. fsolve fails to solve the two equations (really numerical procedures = 0).
But I force it to print out reasonably good solutions.
3. Knowing a reasonably good solution I use that as an initial point for optimizing the sum of squares of the two quantities we want to be zero.
This turns out to give a very good solution.

The worksheet:
MaplePrimes14-10-10odesysAA.mw

Please feel free to ask further questions.