@torabi
I think you just refer to the whole system eq1, eq2, eq3, where the parameter p is present in eq1.
I have no name for the idea of adding an extra condition to determine p. That feature has been available in dsolve for a long time.
Basically, you could do it yourself by replacing the parameter p in eq1 with p(x) and then adding a fourth ode:
eq4:= diff(p(x),x) = 0; # This makes it certain that p(x) will be constant.
So you do:
dsolve({eq1,eq2,eq3,eq4} union bcs union { extra }, numeric);
## But why go the trouble? Maple does that by itself.
 

There are several problems of syntax.
1. You cannot use { } or [ ] as substitutes for parentheses. You must use ( ).
2. It appears that there are missing multiplication signs.
3. By (d(C(t,r))/dt)^2  do you mean (diff( C(t,r), t))^2 (i.e. the first derivative squared) or do you mean the second derivative, i.e. diff( C(t,r),t,t) ?

@9009134 In my answer above I have now at the bottom placed the proof that eq1 and eq2 with bcs23 only has the trivial solution.
The proof just uses the standard procedure (from scratch) for linear bvps with constant coefficients.

In my Maple 17.02 I have no problem with this:

restart;
t0:=time();
int(3628800 / (y * (1 / 2 + y)^11) - 3628800 / (y * (39 / 2 + y)^11), y = 39 / 2..infinity);
time()-t0;

The time is reported as approximately 0.2 sec.
kernelopts(version);
Maple 17.02, X86 64 WINDOWS, Sep 5 2013, Build ID 872941

You should be using fsolve instead of solve if you want to solve numerically.
 

@Tom68 You wrote:

"how to bring to the equation an initialize force? I.e. force of first spring activation?

m*diff(x(t),t,t) + k*x(t) + mu*m*g*signum(v(t)) = F  ???"

To get things started I used the initial conditions x(0) = x0, D(x)(0) = 0, where x0 > 0.
Thus you start at a top in my formulation.
How you accomplish that start is irrelevant to the mathematical question.
 

@vv Thank you. I wonder who wrote the equations shown in the image

It is striking how similar the graphs presented above look like the graphs produced by numerical solution using events as described in my other answer. Notice that params had mu=0.0003, so not as here mu = 0.003. That wasn't intended, but let us stick with the latter one, mu = 0.003.
My conclusion is that the OP's images of the equations ought to have referred to x'(t) being positive or negative, not to x''(t), as I also (luckily) read or misread it.
 

restart;
Digits:=15:
sol:=(A-n*C)*cos(omega*t)+(-1)^n*C/2;
x0:=eval(sol,{n=0,t=0}); #The actual initial value x(0).
x1:=eval(diff(sol,t),{n=0,t=0}); #The initial value of x'(0) is 0.
params:={m=1.0, k=0.1, g=10.0, mu=0.003}:
C:=2*mu*m*g/k;
omega:=sqrt(k/m);
n:=floor(omega*t/Pi);
plot(eval(sol,params union {A=1}),t=0..250,numpoints=10000, size=[800,default],caption="Initial value of A = 1" ); p1:=%:
plot(eval(sol,params union {A=10}),t=0..250, numpoints=10000,size=[800,default],caption="Initial value of A = 10" ); p10:=%:
#### Now the events version:
ode:=m*diff(x(t),t,t)=-k*x(t)-mu*m*g*(2*b(t)-1);
resE1:=dsolve(eval({ode,x(0)=x0,D(x)(0)=0,b(0)=0},params union {A=1}),numeric,
    discrete_variables=[b(t)::boolean],events=[[diff(x(t),t)=0,b(t)=1-b(t)]],abserr=1e-15,relerr=1e-13,maxfun=0);
plots:-odeplot(resE1,[t,x(t)],0..250,size=[800,default],color=blue,numpoints=10000); p1E:=%:
plots:-display(p1,p1E);
resE10:=dsolve(eval({ode,x(0)=x0,D(x)(0)=0,b(0)=0},params union {A=10}),numeric,
    discrete_variables=[b(t)::boolean],events=[[diff(x(t),t)=0,b(t)=1-b(t)]], ,abserr=1e-15,relerr=1e-13,maxfun=0);
plots:-odeplot(resE10,[t,x(t)],0..250,size=[800,default],color=blue,numpoints=10000); p10E:=%:
plots:-display(p10,p10E);

Here I just plant the last image showing the two graphs for A = 10 on top of each other.
There are actually two, but the blue covers the other:

@tomleslie I'm sure that I don't understand the (for me) revised question, where it is the sign of x''(t) that determines the right hand side of the ode.
So we have for x''(t) >0 that

m*diff(x(t),t,t)=-k*x(t)-mu*m*g

thus we have from that equation that -k*x(t) - mu*m*g > 0, i.e.  x(t) < -mu*m*g/k.
Similarly, for x''(t) < 0 we use

m*diff(x(t),t,t)=-k*x(t)+mu*m*g

so in that case x(t) > mu*m*g/k.
Therefore the question: What happens if x(t) is between -mu*m*g/k and mu*m*g/k ? What is the ode?

@tomleslie I misread the tiny images, but after magnifying those low resolution images I see that probably the sign of x'' is determining the sign of the friction rather than x' as I used.
If you want to use events with a trigger containing the second derivative, then you must introduce a second ode since only x and x' are available in numerical solution otherwise.
So something like

restart;
params:=[m=1.0, k=0.1, g=10.0, mu=0.003]:
ode1:=m*diff(x(t),t,t)=-k*x(t)-mu*m*g*(2*b(t)-1);
ode2:=diff(x(t),t,t)=a(t);
resE:=dsolve(eval({ode1,ode2},params) union {x(0)=1,D(x)(0)=0,b(0)=0},numeric,
    discrete_variables=[b(t)::boolean],events=[[a(t)=0,b(t)=1-b(t)]]);
plots:-odeplot(resE,[t,x(t)],0..250,size=[1800,default]);

The problem with this is that oscillations grow.

The reason for using b(0)=0 in my original answer was that this is consistent with a beginning downward movement.

The symbolic solution (with the ode eqn:=m*diff(x(t),t$2)=-k*x(t)-signum(diff(x(t),t$2))*mu*m*g; )
will begin by isolating the second derivative and finds two odes just as in the OP's question. Thus you should not solve with fixed initial conditions because the solutions of the two odes have to be pieced together.

@Scot Gould 

I see all 5 integers when using Maple notation (aka 1D) in input:
 

if true then 1; 2; 3; 4; 5 end if;

I only see the the last integer, i.e. 5, when using 2D math input:

Now, so far this is what we see. Since print is somewhat special I tried assignment:
 

if true then a:=1; b:=2; c:=3; d:=4; e:=5 end if;
a,b,c,d,e;

All works as expected. The 5 assignments are made.
Now 2D math input:

Now the assignments to A, C, and E are made (and printed during execution of the loop), but B and D1 are still unassigned.
Note: Doing the whole thing after a restart I got the same as stated with the exception that only the assignment to E was printed during execution of the loop, but the assignments to A, C, and E were made.

This is horrible!

kernelopts(version);
`Maple 2017.2, X86 64 WINDOWS, Jul 19 2017, Build ID 1247392`

Notes:
1. If true is replaced with 1=1 the 2D math code works.
2. If the code is wrapped in a procedure it works, as in
p := proc () global A, B, C, D1, E; A := 1; B := 2; C := 3; D1 := 4; E := 5; NULL end proc

or using locals:
q := proc () local A, B, C, D1, E; A := 1; B := 2; C := 3; D1 := 4; E := 5; [A, B, C, D1, C] end proc

That is some comfort to those who use 2D math input, which I don't.

@ Yes, for epsilon = 10 we must do something. We can use initmesh=256, maxmesh=2048. It appears to me that HDMadapt does similar things by itself (I just noticed the printouts from HDMadapt in this and other cases).
You wrote:
       "PS, it is not my intention to make this thread a discussion on when dsolve/numeric fails for BVPs and how to fix it. "
I didn't mean to do that either. But since you make claims about the superiorty of your HDM code over that of dsolve/numeric/bvp for certain cases, I found it relevant to ask you to point to a concrete example of that superiority.
You have done that. Thanks.

@ Thank you for the reply.
Using continuation or finding an appropriate approximate solution (different from the one calculated by dsolve/numeric/bvp itself) certainly can require some ingenuity.
In the case you give above with epsilon=8, however, it requires no ingenuity to try increasing maxmesh as that is suggested by the error message from dsolve.  maxmesh = 512 works:

sol1:=dsolve({op(subs(pars,EqODEs)),y1(0)=0,y1(1)=1},numeric,maxmesh=512);

But I shall have a look at the Cash examples you have on your website.

## I had a look and realized that I had taken a close look at the examples earlier in connection with some other testing.
Problem 34 is particularly interesting because there are two solutions for each of the epsilon values considered by Cash including for epsilon=3.5, which you consider. Your procedure finds the smaller one of the two as does dsolve/numeric without a user given approximate solution.
The larger one we can get by doing this (here I have epsilon=3.5):
 

## First the smaller:
resP:=dsolve(eval({EqODEs[],eval(bc1,x=0)[],eval(bc2,x=1)[]},pars),numeric,output=listprocedure);
Y1p,Y2p:=op(subs(resP,[y1(x),y2(x)]));
## Then the larger using an approximate solution:
res2:=dsolve(eval({EqODEs[],eval(bc1,x=0)[],eval(bc2,x=1)[]},pars),numeric,approxsoln=[y1(x)=2*Y1p(x),y2(x)=2*Y2p(x)]);

The gap between the two is larger for smaller values of epsilon.

It could be helpful if you could point to some concrete examples where your code performs better than dsolve/numeric or where the latter fails.

@John SREH It confuses me when I read your statement: " But for other equation the singularity will show up "

A major point in my answer was that the singularity for the given equation always shows up!