@Earl That value of DthetaZero has been found already. It is found above by fsolve and is L[2].
L[1] is the value of t for which this happens. Notice the formal arguments (parameters) for the procedure Q. They are t and DTZ in that order. DTZ was meant to remind us that it refers to DthetaZero.
There is nothing to prevent you from replacing the name DTZ with the name DthetaZero in Q.
The global names t and DTZ (or t and DthetaZero) will never be assigned values by Q.

If you try

fsolve(T,complex,fulldigits);
#or just
fsolve(T,complex);
##  you will see the same behavior. 15 real roots are found by fsolve(T) and fsolve(T,fulldigits).
 

@Thomas Richard After reading your answer I also tested the Maple versions I mention in my reply using Notepad++. The result was the reverse of my tests using a file saved in Notepad (Windows).

Just to include the Danish characters æ, ø, and å I tried this version of your procedure.

test:=proc()
print("hüstel");
print("Rødgrød med fløde på en æblebund");
end proc:
restart;
read("F:/TEST.mpl");
test();
##
The characters ü, æ, ø, and å were replaced by squares in Maple 18 and later.
In Maple 12, 15, 16, and 17 (on the same computer!) the same file worked as it should.
### In the above test the file was saved in Notepad (Windows product).
### When using Notepad++ instead there was no problem in Maple 2017, 2016, 2015, 18, but problems in Maple 12, 15, 16, and 17.

@Christopher2222 I referred to the monstrosity +'ve.
 

@n78745 Looking into kzm and kzp I find that for Vx > 1.116810485*10^6 you should have luck.
I get:

What I did was to execute your whole worksheet (without the plot command obviously and I used Digits:=15 at the top).
After that I did:
 

## To simplify I use that Vx and f only appear as Vx^2 and f^2:
##
kzm2,kzp2:=op(eval([kzm,kzp],{Vx=sqrt(Vx2),f^2=f2})) assuming Vx2>0;
kzm2:=factor(kzm2);
kzp2:=factor(kzp2);
## Looking at kzm2 and kzp2 I notice that qm must be negative and qp positive:
qm:=op([2,1],kzm2);
qm:=eval(numer(qm),f2=1);
qp:=op([2,1],kzp2);
qp:=eval(numer(qp),f2=1);
plot(qm,Vx2=10^11..10^13); 
plot(qp,Vx2=10^11..10^13);
## The minimal value of Vx seems to be:
Vxmin:=sqrt(1.24726565884742*10^12);
## The plot:
plots:-implicitplot(f1,f=10..1e5,Vx=Vxmin..10^8);

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