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!

I don't see any attachment.

I just saved the attached worksheet to the cloud first as public (but that takes a while to show up and probably won't as it is not of general interest) and then as private.
The private one was retrieved and looked like the attached worksheet, i.e. in good shape.
MaplePrimes17-08-23_vokaltest.mw

I realize that I didn't use the workbook format. My test was with a simple worksheet.

So you have a concrete ode of order three, which Maple could easily solve using dsolve:
 

restart;
ode:=diff(y(x),x,x,x)+2*diff(y(x),x,x)-9*diff(y(x),x)-18*y(x)=18*x^2-18*x+22;
dsolve(ode);

Now what is f supposed to represent?
I ask because if you want to write this ode as a first order system (a very common approach), you need to introduce 3 functions f, g, and h (or whatever you would like to call them) with f(x)= y(x), g(x) = diff(y(x),x), and h(x) = diff(y(x),x,x).

Besides, you are talking about an ivp (i.e. an initial value problem), but what are the initial values?
## Solving with initial values using ics as given here:

ics:=y(0)=y0,D(y)(0)=y1,(D@@2)(y)(0)=y2;
dsolve({ode,ics});

I don't think that the idea is viable. Patches are best left to professionals.
We happy amateurs can once in a while find a workaround for some bug or inadequacy. In fact that can be a lot of fun.

Do you envision a version number on, say, patch12? Obviously, the patch idea if successful would make patch12 grow and need modification all the time. By whom? Anybody?
Is your vision a Maple 12 with no errors although with no improvements or features from later versions?

Remember that the update of the Physics package wasn't done by an amateur!
 

You have missing multiplication signs in several places in eq1, eq2, eq3.
When I execute eq1 as you wrote it I get (copying as an image):

If I insert multiplication signs where they are clearly needed I get the input equations:
 

eq1:=(diff(psi(x), x))^2+(diff(u(x), x)+4*(diff(w(x), x))^2)*(diff(psi(x), x))^2+3*(diff(diff(w(x), x), x))+5*(diff(diff(w(x), x), x))*(diff(psi(x), x))-7*(diff(diff(diff(u(x), x), x), x))-28*(diff(diff(w(x), x), x))^2-(21/2)*(diff(diff(diff(w(x), x), x), x))*(diff(w(x), x))+3 = p;
##
eq2:=20*(diff(diff(psi(x), x), x))+50*(diff(diff(diff(w(x), x), x), x))+8*(diff(diff(diff(diff(psi(x), x), x), x), x))-7*(diff(w(x), x))+7*psi(x) = 0;
##
eq3:=4*(diff(diff(w(x), x), x))-4*(diff(psi(x), x))+34*(diff(diff(diff(psi(x), x), x), x))+26*(diff(diff(diff(diff(w(x), x), x), x), x)) = 0;

As you will notice, they are quite different from yours. That is a consequence of the bad decision (in my opinion) of allowing spaces to be interpreted as multiplication in 2D math input. People forget the space occasionally.
I never use 2D math input.
I shall correct my answer below to reflect the revised equations.
I didn't see the problem till today. And why not? Because the original equations also make sense and happen not to become totally ridiculous (that would have saved us).
Example:
In the original eq3 you have a missing multiplication sign in the term (23+11)(diff(psi(x), x, x, x)).
Now this is automatically simplified to 34(diff(psi(x), x, x, x)). In Maple numerical constants like (here) 34 act as constant functions so that 34(y) = 34 for all y no matter what y is. Thus 34(diff(psi(x), x, x, x)) = 34.