@claus christensen I have been using Maple with Windows for many years (since Maple V, Release 2,1992, I believe).
I have had problems with the Danish characters in the past too.

@ThU The question was: how do i select different literal subscripts in a formula in Maple?
The question is not quite clear: Is the meaning: Find the terms with the literal subscript b?
Or does it really have nothing to do with literal subscripts (taken to mean names generated as in xxx__yyy), but is just: Find the first and third operand in the sum.
In the latter case you provide the answer.
But if we understand the question in the first sense, then what to do with something like
u:=a__b+b__a+a__b^2+7+a*b+q__b;

How do you select the terms having xxx__b  for some name xxx without relying on the visual order of the terms?
Since
whattype(a__b);
returns symbol, it is not immediately obvious.

@Kitonum Yes, but what is its advantage then?

@claus christensen If you yourself use different Maple releases at various occasions, or if you are writing a package to be used by other people who might have older Maple releases, then you could have a problem.
I made a test. I wrote two simple .mpl files in Notepad++ (which by default) saves in UTF-8.
The first file TEST.mpl contains the definition of a package (a module) that only exports the procedure to be found in the other file TestProcedure.mpl.
Here are the files:
 

TEST:=module() option package; export TestProcedure;
$include "C:/Users/Bruger/maple/toolbox/TEST/Procedurer/TestProcedure.mpl"
end module;

TestProcedure:=proc()
print("hüstel");
print("Rødgrød med fløde på en æblebund")
end proc:

Then I used my usual way of saving this package as an archive named TEST.mla.
Testing simply like this:

restart;
with(TEST);
TestProcedure();

It worked nicely in Maple 2017.2, but not in Maple 17.02.

@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)?