You are saying that B2211 is positive since B211 and B221 are.
The problem with that is that both B211 and B221 don't remain positive.
MaplePrimes18-03-24_odesys.mw

The plot of B211 and B221:

If your system doesn't behave the way you expect I'm pretty sure it has absolutely nothing to do with odeplot.
My guess is that the problem (if there is one) is in your own formulation of the problem sent to dsolve.
Just take a look at sol(3000) which produces

quite independent of odeplot. I'm also pretty confident that dsolve doesn't make blunders like that.

Nothing is attached.

@CK1992 What it means is that
 

exp(g(x,0)*sqrt(x)/(1+varepsilon*g(x,0)*sqrt(x)))
## is replaced by
1+g(x,0)*sqrt(x)/(1+varepsilon*g(x,0)*sqrt(x)

Starting from BC you can make the replacement like this:
 

BC2:=eval(BC,exp=(x->1+x));

so when varepsilon=0 you get
 

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

I did that in a desperate attempt to get some result. That replacement can only be expected to be any good if the argument to exp is small. The argument is g(x, 0)*sqrt(x)/(1+varepsilon*g(x, 0)*sqrt(x)) and that is only small if g(x, 0)*sqrt(x) is small.
So it doesn't look promising, does it?

@CK1992 In fact I tried using an eta value greater than zero. That worked better at least when using the approximation exp(x) = 1+x. Also I used spacestep=0.1.
Another thing: The boundary value for D(T) is given at eta=20, i.e. you have D(T)(20) = 0.
Could that (or should that) be chosen smaller, e.g. eta=10?
The reason I ask is that the plots

plots:-odeplot(sol1,[[eta,T(eta)],[eta,diff(T(eta),eta)]]);

show that T and T' has become basically constant already at eta=10.
Using that also, this produced a seemingly exponentially growing g*sqrt(x):
 

pds:-plot(g*sqrt(x), x = 0 ..3 , eta = .0001, color = red);

But remember, I was using exp(x) = 1+x, so no great news.

@CK1992 I left out the change in alpha from 1 to zero.
So now a "solution" is produced.
The input to pdsolve when different from the names of the dependent variables (as in this case) should probably include the independent vaiables as in: g(x,eta)*sqrt(x). 
Besides that I didn't have much success. I tried replacing exp by its first taylor polynomial, thus using exp(x) = 1 + x.
Only limited success. 

@CK1992 I tried with these conditions you just gave:
 

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

BCtest := {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}

From
 

pds := pdsolve(PDE1, BCtest, numeric, spacestep = .25);

I get the same error received before: too complex!
No wonder though, because your "new" conditions are exactly the same as the old!
Did you copy the wrong ones?

I have reason to believe that that problem is gone soon.

@Parham2016 Since dsolve uses solve for the job, you might as well be using solve.
One note: Remove the tilde after indets.

indets~(eval({E1, E2},B=10)); #returns {{C1,Nu}}
whereas

indets(eval({E1, E2},B=10)); # returns {C1, Nu}

solve wants the latter form. Apparently the double curlies don't bother dsolve.

Why is there only one solution for Nu, C1 for given B?
If so why must C1 be an increasing function of B?
How can you expect to get dsolve to do anything sensible on {E1,E2} ? They aren't differential equations.
##NOTE: I see that you actually do get a result from dsolve.
It never occurred to me to try something like this:
 

restart;
dsolve({x^2+2,y+2},{x,y});
dsolve({y+2},{y});

The first one works, the second doesn't in Maple 2017.3.
I find this inconsistent. I would like the error in both cases:
Skæg for sig og snot for sig! (I won't translate that).

It appears that you are missing a multiplication sign in the very last statement:
 

(`mod`((`mod`(A, n))*(`mod`(B, n)), n))((`mod`(A, n))*(`mod`(B, n)));

Presumably you wanted:
 

(`mod`((`mod`(A, n))*(`mod`(B, n)), n))*((`mod`(A, n))*(`mod`(B, n)));

Besides that I'm lost.

@adxters If that is what you want then it is much less convoluted:
 

for i to n do 
   if H[i] < 2.7 and A[i, f] < 12 and 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; 
   print(Q[i, foo]) 
end do;

Clearly you will need assumptions. Could you provide text instead of an image and information about r, M, and eta, i.e. are they real, or even positive. Is r>2*M or ... ?

@_Maxim_ Yes indeed. Have you submitted a bug report?

@Adam Ledger Sorry, if my formulation was offensive in any way.
I wouldn't reply as extensively as I have if I didn't find your question very interesting.
As I wrote earlier (see above):
" So the question is rather: why is it good (or even necessary) for _Z to be global in RootOf? "

If _Z had been local, then it would be local to the procedure RootOf (or some other procedure called by RootOf).
That would mean that different calls to RootOf would return different local variants of _Z.
To see that take my example above:
 

restart;
## The home made procedure.
## Remark: Notice that rootof has never been defined.
##
R:=proc(pol,x) local _Z; rootof(subs(x=_Z,pol)) end proc;
##
R(x^2-2,x);
z1:=op(indets(%,name)); # The first local _Z
R(x^2-2,x);
z2:=op(indets(%,name)); # The second local _Z
z1-z2; # NOT zero
## The addresses of the two local variants of _Z are different:
addressof(z1);
addressof(z2);

That fact may be a problem for RootOf.