Just from picking rather arbitrary values for the parameters you get
eval(cdm_ode,{c0=.2, n=5,sigma=175, s0=200, ks=.5, h1=1000, h2=.1,kp=0.5, A=1, B=2});
#It is not very surprising then that
err(.2,5,200,.5,1000,.1,.5,1,2);
detects a singularity.
Trying the presumably simplest case n=1:

sys:=eval(cdm_ode,{c0=.2, n=1,sigma=175, s0=200, ks=.5, h1=1000, h2=.1,kp=0.5, A=1, B=2});
res := dsolve([sys, y1(0) = 0, y2(0) = 0, y3(0) = 0, y4(0) = 0, y5(0) = 0, y6(0) = 175], numeric, range = 0 .. tol_t);
plots:-odeplot(res,[t,y1(t)],0..10);
#The singularity seems obvious.

@Dmitry You cannot know that it always has an additional solution, but you cannot know that it doesn't either.
eq1:=simplify(eval(eq,params=~{seq(1..nops(params))}));
res:=solve(eq1,t3);
evalf[50](eval(res,t1=1));
plots:-implicitplot(eq1,t1=0..5,t3=0..5,gridrefine=3,grid=[50,50]);
plot(eval(lhs(eq1),t1=1),t3=0..2,-1..1);
Digits:=20:
fsolve(eval(eq1,t1=1),t3=1.33);

#So here there are clearly 2 positive solutions.

@Dmitry No, on the contrary. I mean that generally you should expect more than one solution (even more than one real solution). I added a few lines above which shows that when all parameters are 7 and t1=1 then besides t3=1 and t3 = 1.067641139 there is also a complex solution t3 = -1.238855445-.8386023959*I, and even one more real solution.

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]:=...