In your procedure replace nops by numelems.
 

restart;
#Example:
v:= <1, 2, 3, 2, 4>;
nops(v);
numelems(v);
op(v);
for i from 1 to nops(v) do op(i,v) end do;

I shall assume that you only want the answer for sin(x) >=0.
Thus I shall only consider 0<= T <= Pi.
 

restart;
res:=int(sqrt(sin(x)),x=0..T) assuming T>0, T<Pi;
diff(res,T);
simplify(%) assuming T>0,T<Pi; #WRONG on Pi/2..Pi
## Thus we split:
res1:=int(sqrt(sin(x)),x=0..T) assuming T>0, T<Pi/2;
simplify(diff(res1,T)) assuming T>0,T<Pi/2; #OK
res2:=int(sqrt(sin(x)),x=Pi/2..T) assuming T>Pi/2,T<Pi;
simplify(diff(res2,T)) assuming T>Pi/2,T<Pi; #OK
resp:=piecewise(T<Pi/2,res1,eval(res1,T=Pi/2)+res2); # Result on 0..Pi
plot(resp,T=0..Pi);

The title of your question was: "More than one solution for a special case."
But fsolve will only give you one result unless the input is a polynomial in one variable.
###
Well, since dsolve gives you a result, try it again with a higher setting of Digits.
Just for a short experiment I set Digits to 30 (you can try much higher):
Removing the imaginary solutions I got these 4 solutions (when B=10):
[{C1 = -8.21367014680801194633270213292*10^(-12), Nu = -4.38002192864507929874554229402*10^6}, {C1 = 20.4996314736857070747869876646, Nu = -65.8211157946078643801270593845}, {C1 = -.898697538127575625586802649353, Nu = -53.0523912091922304391249719870}, {C1 = -8.21426028148815833335058821461*10^(-12), Nu = 4.37969614250226629250101130330*10^6}].
###########
Trying your loop using dsolve instead of fsolve and Digits=50 and also removing imaginary results.
You might as well be using solve. The use of dsolve here is confusing (was to me!).

Bvals:=[ -10.0, -1.0, -0.5, 0, 0.5, 1.0, 10.0
       ]:
for j from 1 to nops(Bvals) do
   temp:= [evalf(dsolve(eval({E1, E2},B=Bvals[j]), indets(eval({E1, E2},B=Bvals[j]))))];
   res[j]:=remove(has,temp,I)
end do:
for j from 1 to nops(Bvals) do nops(res[j]) end do;

I'm getting 2 real results for the first 4 B values and 4 for the last 3 (edited).
I would suggest using solve instead of dsolve.

Isn't this just a case of garbage in, garbage out?
Using allvalues gives me 2 matrices:
 

restart;
with(LinearAlgebra):
x := RootOf(z^2-t, z);
m := Matrix(2, (i,j) -> add(a[j, k]*((-1)^(i-1))^k*x^k, k = 0..1));
m1,m2:=allvalues(m);
L1,V1 := Eigenvectors(m1);
## and
L2,V2 := Eigenvectors(m2);

The two outer ifs have no else, therefore if not H[i] < 2.7 or if not A[i, f] < 12 then nothing is going to happen.

So you must tell us (or yourself) what you want to do in those cases.
For clarity I wrote the code with some indentation:
 

for i to n do 
   if H[i] < 2.7 then 
     if A[i, f] < 12 then 
       if A[i, o] < 1.2 then 
         Q[i, foo] := evalf(610*(A[i, o]*sqrt(H[i])*h[k]*A[i, T])^(1/2)) 
       else 
         Q[i, foo] := evalf(7.8*A[i, t]+378*A[i, o]*sqrt(H[i])) 
       end if 
     end if 
   end if; 
   print(Q[i, foo]) 
end do;

You seem to have some syntactical problems with your derivatives. Being no friend of 2D math input I made use of Kitonum's translation into the trustworthy 1D math input, known also as Maple Input.
 

restart;
Digits:=15:
interface(imaginaryunit = j);
k := .5; tau := .95; mu := 0.1e-1; pi := 116.1; theta := 0.8e-2; phi := 0.25e-2; epsilon := 0.2e-2; rho := 0.5e-1; beta := 0.115e-1; chi := 0.598e-2; q := .5; eta := .2; delta := .1; alpha := 0.57e-1; p := .2; Upsilon := 1.2;
lambda := k*tau*(C(t)*Upsilon+(I)(t))/(S(t)+V(t)+C(t)+(I)(t)+R(t));
sys:={diff(C(t), t) = rho*lambda*S(t)+rho*epsilon*lambda*V(t)+(1-q)*eta*I(t)-(mu+beta+chi)*C(t), diff(I(t), t) = (1-rho)*lambda*S(t)+(1-rho)*epsilon*lambda*V(t)+chi*C(t)-(mu+alpha+eta)*I(t), diff(R(t), t) = beta*C(t)+q*eta*I(t)-(mu+delta)*R(t), diff(S(t), t) = (1-p)*pi+phi*V(t)+delta*R(t)-(mu+lambda+theta)*S(t), diff(V(t), t) = p*pi+theta*S(t)-(epsilon*lambda+mu+phi)*V(t)};
ics:= {S(0) = 8200, V(0) = 2800, C(0) = 200, (I)(0) = 210, R(0) = 200};
res:=dsolve(sys union ics, numeric); 
plots:-odeplot(res,[C,I,S](t),2000..3000,numpoints=10000);
plots:-display(Array(1..5,1..1,[seq(plots:-odeplot(res,[t,F(t)],2000..3000,numpoints=10000),F=[C,I,R,S,V])]));

The 3 dimensional plot of C,I,S is shown here. Notice that I only plot after the initial settling has taken place. I chose 2000..3000:

Below we see that your system has 3 equilibrium points and that they all are unstable.
 

RHS:=convert(evalindets(rhs~(sys),function,x->op(0,x)),list);
sol:=[solve(RHS,{C,I,R,S,V})];
pts:=map(subs,sol,[C,I,R,S,V]);
J:=unapply(VectorCalculus:-Jacobian(RHS,[C,I,R,S,V]),C,I,R,S,V);
for pp in pts do simplify(LinearAlgebra:-Eigenvalues(J(op(pp)))) end do;

In maple z^(1/3) is the principal cube root of z. Thus as an example (-1)^(1/3) is equal to 1/2+I*sqrt(3)/2.
Try evalc( (-1)^(1/3) );
What you can use since you expect to get -1 from (-1)^(1/3) is surd.

Thus your plot will be produced by
 

plot(surd(x^3-x^2,3),x=-3..3);

Your code has some syntactical problems. I tried to clear them up.
Besides that the second derivatives w.r.t. eta cannot appear in the boundary conditions. I removed them.
Thus the code looks like this:
 

restart;
Digits := 10;
sys_ode := 2*diff(T(eta), eta, eta, eta)+T(eta)*diff(T(eta), eta, eta) = 0;
ics := T(0) = 0, D(T)(0) = 1, D(T)(20) = 0:
sol1 := dsolve({ics, sys_ode}, numeric, output = operator);
q,w:=op(subs(sol1,[T,D(T)]));
######
PDE1 := eval([2/P * diff(g(x, eta), eta, eta)+q(eta)*diff(g(x, eta), eta)-g(x, eta)*w(eta) = 2*x*w(eta)*diff(g(x, eta), x), 2/S * diff(phi(x, eta), eta, eta)+q(eta)*diff(phi(x, eta), eta) = 2*x*w(eta)*diff(phi(x, eta), x)]);
##
subBC1 := -phi(x, 0)*exp(g(x,0)*sqrt(x)/(1+varepsilon*g(x,0)*sqrt(x)));
subBC2 := alpha*phi(x, 0)*sqrt(x)*exp(g(x,0)*sqrt(x)/(1+varepsilon*g(x,0)*sqrt(x)));
##
BC := {g(0, eta) = 0, g(x, 20) = 0, phi(0, eta) = 1, phi(x, 20) = 1, D[2](g)(x, 0) = subBC1, D[2](phi)(x, 0) = subBC2};
##
P := 1;
S := 1;
varepsilon:=1: #Added
alpha:=1: #Added
#######
pds := pdsolve(PDE1, BC, numeric, spacestep = .25);

I receive the error:

Error, (in pdsolve/numeric/match_PDEs_BCs) unable to associate boundary conditions and dependent variables, system is too complex
Even when setting varepsilon=0 the same error comes up.
In an test to see if there may be other problems I tried these boundary conditions, where exp is replaced by 1:
 

BCtest:={g(0, eta) = 0, g(x, 20) = 0, phi(0, eta) = 1, phi(x, 20) = 1, D[2](g)(x, 0) = -phi(x, 0), D[2](phi)(x, 0) = alpha*sqrt(x)*phi(x, 0)};

The answer came up right away with no complaint, so it appears that the complexity of BC is the problem.

You could try simplify/size:
 

length(Teller); # 4308
TS:=simplify(Teller,size);
length(TS); # 2035
whattype(TS); # `+`
nops(TS); # 5

Begin by not assigning to the parameters you want to change later.
 

restart;
lambda := 1.0; K__r := 1.0; S__r := .5; m := .5; M := sqrt(10.0); `ϰ` := .5; Omega := sqrt(5.0); 
## Have omitted Gm, Gr, P__r, S__c
## Now define PDE and IBC as you did
##  ...
## Define SOL as a procedure:
SOL:=(gm,gr,pr,sc)->pdsolve(eval(PDE,{Gm=gm,Gr=gr,P__r=pr,S__c=sc}), IBC, numeric, spacestep = 0.1e-1);
## Test:
SOL(5,6,0.71,0.22):-plot(t = .3, color = red);
SOL(5,6,0.71,0.1):-plot(t = 1,color = red);

Adding this code you can get an animation;
 

PL1:=sc->SOL(5,6,0.71,sc):-plot(t = 1,color = red);
animate(PL1,[sc],sc=0.01..1);

I'm assuming that you meant sys:= {...} and not sys*{...}.
Thus we have:
 

sys:={A2+D2 = 0, B1*sin(192*K1)+D1 = 0, 3.383*B2*K1+C2 = 0, B1*K1^2*sin(192*K1)+11.444689*A2*K1^2*cos(568.344*K1)+11.444689*B2*K1^2*sin(568.344*K1) = 0, A2*cos(568.344*K1)+B2*sin(568.344*K1)+168*C2+D2 = 0, -3.383*A2*K1*sin(568.344*K1)+3.383*B2*K1*cos(568.344*K1)+C2-B1*K1*cos(192*K1) = 0};
##You cannot use fsolve if you want to keep K1 unassigned.
solve(sys,{A2, B1, B2, C2, D1, D2});

You just get the trivial solution {A2 = 0., B1 = 0., B2 = 0., C2 = 0., D1 = 0., D2 = 0.}.
## You may ask: For which values of K1 are there nontrivial solutions?
Notice that your equations are linear in A2, B1, B2, C2, D1, D2.
To look into that use a little linear algebra:
 

A,zz:=LinearAlgebra:-GenerateMatrix(sys,[A2, B1, B2, C2, D1, D2]);
dt:=LinearAlgebra:-Determinant(A); #Must be zero to get nontrivial solutions
plot(dt,K1=-.1..0.1);
plot(dt,K1=-.01..0.01);
fsolve(dt=0,K1=0); # Answer 0.
fsolve(dt=0,K1=0.006); # Answer: 0.5941860004e-2

We must expect infinitely many values of K1 with that property.
# Since dt-eval(dt,K1=-K1) evaluates to zero, you need only look at K1 >= 0.
To find the corresponding solutions for A2, B1, B2, C2, D1, D2 you could do:

K1a:=fsolve(dt=0,K1=0.006); #Example only
G:=LinearAlgebra:-GaussianElimination(eval(A,K1=K1a));
G0:=fnormal~(G); #When dt=0 the rank < full. Removing roundoff.
LinearAlgebra:-LinearSolve(G0,zz,free=t);

S:=cos(440*2*Pi*t)+cos(554*2*Pi*t)+cos(659*2*Pi*t);
convert(S,exp);

I changed pi to Pi.

This probably won't be satisfactory, but I tried without using plotsetup. I right clicked on the plot and exported as a jpg file:
 

plots:-logplot(x^2-3, x = 0 .. 100, font = [bold, "TimesNewRoman", 30],size=[1500,1100]);

The plot:

The functions named _F1, _F2, ... are meant as arbitrary functions. Thus an expression like _F1(t+x) means any expression depending on the sum t+x only. Nothing more is meant (except that appropriate smoothness requirements are implicitly assumed).

You need initial conditions for the derivative too. Your ode is called ode1, not as you have ode.
To get arrows you must rewrite as a system of first order:
 

restart;
ode1 := diff(x(t), t, t) = (5*9.80665)*sin((1/6)*Pi)-(10*(10-sqrt(x(t)^2+25)))*x(t)/sqrt(x(t)^2+25)-(diff(x(t), t));
ics:=[seq([x(0) = (1/4)*k,D(x)(0)=k/2], k = -20 .. 20)]; # I chose D(x)(0)=k/2
DEtools[DEplot](ode1, x(t), t = -2 .. 2, ics, x=-7..20, stepsize = 0.5e-1, linecolour = red);
ode2:=diff(x(t),t)=y(t);
sys := {ode2,subs(ode2,ode1)};
ics_xy:=subs(D(x)=y,ics);
DEtools[DEplot](sys, [x(t),y(t)], t = -2 .. 2, ics_xy, x=-7..15,y=-7..20,color = blue, stepsize = 0.5e-1, linecolour = red,arrows=medium);

The last plot: