@awass As Carl already has said, this is the change: local x:=evalf(xx) relative to his.
Now why do that?
Because you have the call

plots:-odeplot(nans, [[y, theta]], y = 0 .. 6);

where

theta:=at(f(y)-(1/2)*Pi, diff(f(y), y));

Thus the first argument to at is f(y)-(1/2)*Pi. That would cause problems in an if statement. Try this e.g.
 

if Pi>0 then print("Pi is greater than zero") else print("xxxx") end if;

The simplest solution to this is the use of evalf. In the present context of numerical solution to odes where in fact the argument f(y) - Pi/2 will be replaced by a sequence of numerical values for f(y) in at when handled by odeplot I see no reason for doing anything else. Some people may suggest using is, as in

if is(0 < x) then arctan(yy/x) elif is(x < 0) then Pi+arctan(yy/x) elif is(x = 0) then (1/2)*Pi end if

(and not using evalf), but I wouldn't.
 

### NOTE about floating points and Pi.
In the most recent versions floating point contagion affects Pi.
In Maple 2017 and 2016 this is the default.
Thus if you haven't set kernelopts(floatPi=false) then Pi + 1.2345 would evaluate directly to a float:
 

kernelopts(floatPi); 
# The output will be true if you are using Maple 2016 or 2017 and have the default setting.
kernelopts(floatPi=false); # This will set it to false (but will return what it was before).
## Now try:
Pi/2+1.2345;

I have the following statement in my maple.ini file:
 

try kernelopts(floatPi=false) catch: end try:

This means that in all Maple versions Pi/2+1.2345 will return Pi/2+1.2345 as I prefer.
The reason for the try statement is that the floatPi option doesn't exist in Maple versions 18 and earlier so the statement kernelopts(floatPi=false) will result in an error. That error is caught by the try statement.
A final note: There is no setting to handle other constants or roots like sqrt(2), so the special handling of Pi I find strange. Try this with either setting of kernelopts(floatPi):
sqrt(2)+1.2345;
sqrt(2)*1.2345;
results sqrt(2)+1.2345 and sqrt(2)*1.2345.

@Markiyan Hirnyk Yes Markiyan, you got a point. The result ought to be correct in the interval x= -Pi/3..Pi/3.
But it is not.
I tried using IntegrationTools:-Change and got a different result, but one which indeed is correct in the interval -Pi/3..Pi/3.
 

restart;
f:=(1+cos(3*x))^(3/2);
A:=Int(f,x);
IntegrationTools:-Change(A,u=3*x/2);
IntegrationTools:-Change(%,u1=sin(u));
expand(%);
simplify(%) assuming u1^2<1;
value(%);
eval(%,u1=sin(u));
F1:=eval(%,u=3*x/2);
plot(diff(F1,x)-f,x=-Pi/3..Pi/3); # zero

The result from IntTutor is
F := (4/9)*sqrt(2)*sin((3/2)*x)^3-(4/3)*sqrt(2)*sin((3/2)*x);
That simply has the opposite sign of F1 as found above!

Student:-Calculus1:-IntTutor is obviously designed for students at the level of Calculus 1.
Thus it is a learning tool, not really meant for computational purposes.
The steps taken must be carefully considered by the student and should be considered his responsibility. It is part of his learning process.
Here the step that requires some thinking or caution is the change of variable step u1 = sin(u), where u is the old variable and u1 is the new.  Since sin is not monotone on R, but only piecewise monotone, we cannot expect the final answer to be correct on all of R after getting back to u (i.e. what is referred to as 'revert').
In fact the final answer found is correct on infinitely many intervals, which is illustrated by your plot command but with the range changed:

plot(diff((4/9)*sqrt(2)*sin((3/2)*x)^3-(4/3)*sqrt(2)*sin((3/2)*x), x)-(1+cos(3*x))^(3/2), x = 0 ..4*Pi);

A simpler example to illustrate my point:

restart;
A:=Int(cos(2*x),x);
IntegrationTools:-Change(A,t=sin(x));
## Notice arcsin appearing in the result: Caution needed.
## arcsin is the inverse of sin restricted to [-Pi/2,Pi/2].
expand(%);
simplify(%) assuming t^2<1;
value(%);
F:=eval(%,t=sin(x));
simplify(%) assuming real; # Correct on intervals where cos(x) >= 0 only.

For an example where the answer is only correct on [-Pi/2..Pi/2] try
A:=Int(2*x,x);
with the same change of variable t = sin(x).

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