@9009134 I have no problem executing your statements till and including the ones involving PolynomialInterPolation.
I left out the solve statements since their results are not assigned to, and not used. When I had similar solve statements in the link given above, they where there for the purpose of seeing if solution for the highest derivatives was possible. That was not the case before differentiation, but was the case after.
However, what follows after the PolynomialInterpolation must be missing something or is left over from your previous attempts. There I see a reference to dsys3, which is never defined in that worksheet.

The plots
plots:-odeplot(res,[[x,f1(x)],[x,f2(x)],[x,f3(x)]],0..1);
and
plot([fy11,fy22,fy33],x=0..1);

are (not surprisingly) very similar (indistinguishable). I suppose that your reason for doing the interpolation is that you want analytical expressions for the approximate solutions.

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