I expect that Paul Halmos would have a fit if he saw +'ve xyz axis and wouldn't ever get to the next occurrence of that monstrosity.

@saba Two comments:
1. When you say that " Z(r) is function of integration" do you mean that is an arbitrary function of r (but not of t)?
By that I mean is it similar to the Z(r) I have added as an arbitrary function of r in an attempt to get all solutions of eq:
 

eq:=diff(u(t,r),t)=r*sin(t);
map(int,eq,t) + (0=Z(r));

2. I'm pretty sure that your pde cannot be solved analytically. If I'm right then you would have to solve it numerically.
But that seems out of the question in your case.
To illustrate the difficulties just consider the ode version I gave. Furthermore pick extremely simple concrete values for the quantities independent of t, viz. Z(r), axi, alpha, and c. I picked these:
Z(r)=2,axi=0,alpha=1,c=0.
Then try dsolve on the resulting ode:

restart:
EQ:=diff(C(t,r),t$2)=(Z(r)^2-1)-(C(t,r)^2)/(4*alpha)*(1+((1+(16*(alpha)^2/C(t,r)^4)*(Z(r)^2-1)^2-(4*alpha*c/3)+(8*alpha*axi/C(t,r)^4)))^(1/2));
ode:=eval(EQ,C(t,r)=C(t));
## That was as done in an earlier reply.
## Now use the very simple values mentioned above to get a simpler ode:
ode0:=eval(ode,{Z(r)=2,axi=0,alpha=1,c=0});
## Try to solve that:
dsolve(ode0);
## You are probably not very happy with that result!
## Solve this concrete ode numerically. You need initial (or boundary) conditions. I pick them arbitrarily just for illustration:
res:=dsolve({ode0,C(0)=1,D(C)(0)=3},numeric);
## Plot the graph of C:
plots:-odeplot(res,[t,C(t)],0..5);

@Carl Love Although I have known about coerce I haven't been using it. I should start doing that. Very nice.

I deleted your repost of the question. In that you had (and I quote):

@saba What is Z(r)? Is that a known function?
Notice that your pde is at it stands right now just an ode of second order in t containing a parameter r.
If mehdi jafari's interpretation of your pde is correct then you could turn the pde into an ode by doing:
 

restart:
EQ:=diff(C(t,r),t$2)=(Z(r)^2-1)-(C(t,r)^2)/(4*alpha)*(1+((1+(16*(alpha)^2/C(t,r)^4)*(Z(r)^2-1)^2-(4*alpha*c/3)+(8*alpha*axi/C(t,r)^4)))^(1/2));
ode:=eval(EQ,C(t,r)=C(t));
dsolve(ode);

But surely there must be more to it. What about that Z(r)?

@Markiyan Hirnyk I didn't claim that F1 as I found it above is an antiderivative to f:=(1+cos(3*x))^(3/2) on all of R.
I claimed it was on [-Pi/3, Pi/3]. That it also happens to be an antiderivative on infinitely many other intervals is just due to the periodicity of cos.
I don't think that the ambition of IntTutor is to provide an antiderivative on all of R. At the Calculus 1 level that might not be a good idea.
But at least the various steps performed should be correct. The error in this case lies in the output from IntTutor after the change of variable u1 = sin(u), where the result is:
2*(Int(2*sqrt(2)*u1^2-2*sqrt(2), u1))*(1/3)
The sign of this is simply wrong. It doesn't appear to be due to the use of a different branch of arcsin (like the inverse of sin restricted to [Pi/2, 3*Pi/2] where sin is decreasing). In fact I tried IntTutor again one step at a time, but this time using the change of variable u = arcsin(u1). The output from the change step was the exact same as when using u1 = sin(u).

 

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