@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.

@torabi I chose the condition (D@@2)(u)(0) = 0 as the first simple one I could think about.
The quantities that can be involved in the boundary conditions for your system are all derivatives of orders from j = zero to one less than the highest appearing for each function.
Thus the quanties that can be involved in boundary conditions for your system are:
 

N:=[4,3,4]:
nms:=[psi,u,w];
pts:=[0,1];
S:={seq(seq(seq((D@@i)(nms[j])(a),i=0..N[j]-1),j=1..nops(N)),a=pts)};

Now you already have 11 boundary conditions and they are all of the simple form (D@@i)(f)(a) = c for some i ,f, and c.
So if you want to add an extra simple condition of the same form you can do:
 

## The bcs you have given we shall name bcs:
bcs:={psi(0) = 0, psi(1) = 0, u(0) = 0, u(1) = 0, w(0) = 0, w(1) = 0, D(u)(1) = 0, D(psi)(0) = 0, D(psi)(1) = 0, D(w)(0) = 0, D(w)(1) = 0};
## So you can choose from
S minus lhs~(bcs);

Pick one of those and equate it to some concrete value c, add it to bcs, and use dsolve.
Be aware that p will depend on what you end up doing!
########## To illustrate the last remark try all of this:
 

S1:=S minus lhs~(bcs);
for s1 in S1 do
  try
    res[s1]:=dsolve({eq1,eq2,eq3} union bcs union {s1=0},numeric)
  catch:
    printf(cat(StringTools:-FormatMessage(lastexception[2..-1]),"\n"))
  end try
end do;  
indices(res); # The successes
## We have revised the p results to reflect the revised eqs:
## In all 3 cases we just get the trivial solution
res[(D@@2)(u)(0)](0); # p = 3
res[(D@@2)(u)(1)](0); # p = 3
res[D(u)(0)](0); # p = 3
## Now set the quantities in S1 equal to 1 instead of 0:
for s1 in S1 do
  try
    res1[s1]:=dsolve({eq1,eq2,eq3} union bcs union {s1=1},numeric)
  catch:
    printf(cat(StringTools:-FormatMessage(lastexception[2..-1]),"\n"))
  end try
end do;
indices(res1); # The successes (the same as before)
res1[(D@@2)(u)(0)](0); # p = 13.5000000390248 (the revised eqs)
res1[(D@@2)(u)(1)](0); # p = -17.9999999849048 (rev. eqs)
res1[D(u)(0)](0); # p = -39.0000000347809 (rev. eqs)

If you do the same loop with s1 = 2 you get the following p-values: p = 24, p = -39, p = -81.
So you should ask yourself what you are trying to do!
## If you want to see the graphs of e.g. u in the different cases you can do:

plots:-display(Array([seq(plots:-odeplot(res[rr],[x,u(x)]),rr=indices(res,nolist))]));
plots:-display(Array([seq(plots:-odeplot(res1[rr],[x,u(x)]),rr=indices(res1,nolist))]));

@torabi I don't know what you mean by " is it possible for solving boundary condition.? ".
You haven't commented on my statement " If your intention is to determine p at the same time you solve ...."
I take that to mean that you do want to determine p from the boundary conditions. But then you must come up with an additional boundary condition.
As an example I added the condition (D@@2)(u)(0)=0 besides the one you just gave above: D(u)(1)=0:
 

dsys3a := {eq2, eq3, eq1, psi(0) = 0, psi(1) = 0, u(0) = 0, u(1) = 0, w(0) = 0, w(1) = 0, D(u)(1) = 0, D(psi)(0) = 0, D(psi)(1) = 0, D(w)(0) = 0, D(w)(1) = 0,(D@@2)(u)(0)=0};
res:=dsolve(dsys3a,numeric);
res(0); # you will see that  p = 3.00000000000000. #Using the revised equations: See above
plots:-odeplot(res,[x,u(x)],0..1); #Notice that u is very, very small.

The solution for u, psi, and w above is really the trivial one u = w= psi = 0 identically.
If on the other hand you don't want to impose other boundary conditions than the 11 you got, then you must supply the value of p as I have done in these lines where  I took p = 1:
 

dsys3b := {eq2, eq3, eval(eq1, p = 1), psi(0) = 0, psi(1) = 0, u(0) = 0, u(1) = 0, w(0) = 0, w(1) = 0, D(u)(1) = 0, D(psi)(0) = 0, D(psi)(1) = 0, D(w)(0) = 0, D(w)(1) = 0};
res1:=dsolve(dsys3b,numeric);
plots:-odeplot(res1,[x,u(x)],0..1);