This is tricky. That is why Physics:-Assume exists, I believe.
 

restart;
assume(A < B);
S := 2/(B-A);              
about(A);
about(B);
A;
subs(A=5,B=10,S); # 2/5
################
restart;
Physics:-Assume(A < B);
S := 2/(B-A); 
about(A);
A:=5; B:=10;
S; # 2/5

In view of the result of solving for the derivatives, it is a weird warning:
 

restart;
ode:={diff(x__1(t),t)*sin(x__2(t))=x__4(t)*sin(x__3(t))+x__5(t)*cos(x__3(t)),diff(x__2(t),t)= x__4(t)*cos(x__3(t))-x__5(t)*sin(x__3(t)),diff(x__3(t),t)+diff(x__1(t),t)*cos(x__2(t))= 1,diff(x__4(t),t)-(1-B)*x__5(t)= sin(x__2(t))*cos(x__3(t)),diff(x__5(t),t)+(1-B)*a*x__4(t)=sin(x__2(t))*sin(x__3(t))};
sys:=solve(ode,{seq(diff(cat(x__,i)(t),t),i=1..5)});
dsolve(sys); # No warning (and likely, no result).

Numerical solution causes no problem.
Note:  Setting infolevel[dsolve]:=5: reveals some of what is going on.

Instead of using surd, you could try RealDomain:

with(RealDomain);
plot(x^(1/3),x=-2..2,thickness=3);
plot(x^(1/3),x=-2..2,thickness=3,adaptive=true);

I did this in Maple 2022.0 and got these two graphs:

and

In Maple 2022 adaptive=true isn't the default. In this case it appears that adaptive=geometric (new in 2022) is used for the default.
It doesn't work so well in this case.
Even surd doesn't work well if you stay within RealDomain.
But if you leave RealDomain it works fine:
 

restart;
with(RealDomain);
plot(x^(1/3),x=-2..2,thickness=3); # Not good
plot(x^(1/3),x=-2..2,thickness=3,adaptive=true); # OK
plot(x^(1/3),x=-2..2,thickness=3,adaptive=geometric); # Not good
plot(surd(x,3),x=-2..2,thickness=3); # Not good
unwith(RealDomain);
plot(surd(x,3),x=-2..2,thickness=3); #OK

I shall report this, although it might already be known.

Notice that the help page for SolveEquations says:
"The SolveEquations command can return not only exact solutions of the equation system but also any minimums of function F. If the residuals are too large then the solution is not exact solution, it is rather the solution that minimizes the residuals."

Thus you should throw results with large residuals out:
 

Digits:=15:
A:=SolveEquations(eq, AllSolutions);

select(has,convert(fnormal(A[..,1]),list),0.);
numelems(%);  # In my recent run in Maple 2022.0 I get  18.

So rather than finding too many roots it finds too few. Graphically we see 22 roots.

Also we find 22 roots here:
 

rts:={RootFinding:-Analytic(eq,x=-35-0.1*I..35+0.1*I)};
numelems(rts);

And also 22 here:
 

A:=SolveEquations(eq,[x=-35..35], AllSolutions);

select(has,convert(fnormal(A[..,1]),list),0.);
numelems(%); ## 22

rk2 is described in the help page for dsolve/numeric/classical.
Notice that f is a vector function. Also k1, k2, and Y[n] are vectors in the systems case.
First it solves for the derivatives, then it uses the algorithm described in the help page.
Your system after solving for the derivatives is extremely simple:
 

restart;
odes1 := {diff(A[0](t), t) + 1/3*diff(A[1](t), t) + 1/9*diff(A[2](t), t) - 2*A[2](t) = 0, diff(A[0](t), t) + 2/3*diff(A[1](t), t) + 4/9*diff(A[2](t), t) - 2*A[2](t) = 0, diff(A[0](t), t) + diff(A[1](t), t) + diff(A[2](t), t) - 2*A[2](t) = 0, A[0](0) = -9, A[1](0) = -5, A[2](0) = 5};
odes,ics:=selectremove(has,odes1,diff);
sys:=solve(odes,diff~({A[0],A[1],A[2]}(t),t)); # Very simple
dsolve(odes1);            # exact sol
dsolve(sys union ics);    # exact sol
interface(rtablesize=20);
res:=dsolve(sys union ics,numeric,method=classical[rk2],
            stepsize=0.1,output=Array([seq(k/10,k=0..10)]));
dsolve(odes1,numeric,method=classical[rk2],
            stepsize=0.1,output=Array([seq(k/10,k=0..10)]));

Your ode has 2 equilbria E1 and E2. It appears that the solutions approach those in finite time. That is possible because of the lack of uniqueness at the initial values E1 and E2. This can cause numerical problems when those values are approached,  as it does in this case.

restart;

Digits:=15:
V:=n^2*((T[e]+T[i])/T[e])*((((n-1)*(1+1/(2*delta[p]))))-ln(n)-ln(n/(2*delta[p])));
W:=eval(V,{T[e]=6/5,T[i]=1/100,delta[p]=1/5,n=n(x)}); # Preferring exact values to start with
ode:=diff(n(x),x)^2+2*W=0;
ode1,ode2:=solve(ode,{diff(n(x),x)}); # Two separate odes
E:=solve(rhs~(ode1),n(x));
E1:=rhs(op(E[1])); # Equilibrium 1
E2:=rhs(op(E[2])); # Equilibrium 2 
solve(rhs~(ode2),n(x));# Same values E1 and E2
evalf([E1,E2]);
ic:=n(0)=1;
res1:=dsolve(ode1 union {ic},numeric);
plots:-odeplot(res1,[x,n(x)],-1.9..0.9);
res1(0.9);
evalf(E2);
res1(-1.8085215);
evalf(E1);
res2:=dsolve(ode2 union {ic},numeric);
plots:-odeplot(res2,[x,n(x)],-0.52455902..1.8085215);
res2(1.8085215);
evalf(E1);
res2(-0.52455902);
evalf(E2);
res1(0.524559);

Your equation eq2 is
 

C1*C2*C3*C4(C1*C2*R3*R4 + C1*C3*R2*R4 + C1*C4*R2*R3 + C2*C3*R1*R4 + C2*C4*R1*R3 + C3*C4*R1*R2) = 126

but should be
 

C1*C2*C3*C4*(C1*C2*R3*R4 + C1*C3*R2*R4 + C1*C4*R2*R3 + C2*C3*R1*R4 + C2*C4*R1*R3 + C3*C4*R1*R2) = 126

there may be more `*´ missing in other equations. I strongly advocate using 1D math input. That requires explicit use of `* ` and it is much easier to detect errors like that.

I'm not saying that fsolve will give a solution after that problem is resolved, but it certainly must be resolved first.
##### Added:
You can find where `*` is missing by doing:
 

indets~([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9],function);

The output is:
[{C4(C1*R2*R3*R4 + C2*R1*R3*R4 + C3*R1*R2*R4 + C4*R1*R2*R3)}, {C4(C1*C2*R3*R4 + C1*C3*R2*R4 + C1*C4*R2*R3 + C2*C3*R1*R4 + C2*C4*R1*R3 + C3*C4*R1*R2)}, {}, {}, {}, {}, {}, {}, {}].
This means that the problem exists in eq1 and eq2, but not in the rest.
C4(....) is considered a function, whereas C4*(...) is a product.

Add this as the last argument to polarplot:
coordinateview = [0 .. 1, default]

You could do either one of these instead.  (a third and better way at bottom).
 

restart;
eq:=P=A+(z+x/2/z);
(rhs=lhs)(-(isolate(eq,A)-P));
restart;
eq2:=P=A+``(z+x/2/z); # Notice `` 
isolate(eq2,``(z+x/2/z));
expand(%);

PS. Notice that the help page for isolate indeed does say:

"The procedure isolate attempts to isolate the second argument expr in the first argument eqn and solves eqn for expr."
here expr is "any algebraic expression".
That is rather vague.
But look at the examples on the help page. In all cases expr is an actual operand of one of the sides of the equation.
In your case (z+x/2/z) is not an operand of P=A+(z+x/2/z).
Notice the output when you define the equation: The panthesis is gone:
 

restart;
eq:=P=A+(z+x/2/z);
op(rhs(eq));

Answer: A, z, 1/2*x/z.
Thus a third way around this would be:
 

isolate(eq,x/z/2) + z;

I made a few changes. The main one is that I don't in fact ever assign to S.
In any case you should wait with that until the symbolic computations are over.

I make 'solution' into a procedure, but that is just for ease of comparison of results for the exact 1 and the float 1., which is the point of all this.

restart;
Digits:=15:
#S:=1; # Don't assign to S
F:=add((q^i)*f[i](eta),i=0..16);

eqa:=simplify(diff(F,eta,eta,eta)+(q*g*((F*diff(F,eta,eta))-((diff(F,eta))*(diff(F,eta)))-(S*diff(F,eta))-(((S*eta)/2)*diff(F,eta,eta))))):
f[0](eta):=(alpha/2)*eta^2+eta:

for m from 1 to 8 do eq[m]:=coeff(eqa,q^m) end do:
for m from 1 to 8 do (simplify(int(int(int(((-1)*eq[m]),eta),eta),eta)+f[m](eta))):f[m](eta):=simplify(%) end do:

a1:=simplify(add(f[n](eta),n=0..8)):
b1:=simplify(subs(eta=1,a1))=(S/2):b2:=simplify(subs(eta=1,diff(a1,eta,eta)))=0:
sys:={b1,b2}:

solution:=s->fsolve(eval(sys,S=s));

AG_vals:=[solution(1),solution(1.)]; # (Almost) exactly the same
c1E:=eval(a1,AG_vals[1] union {S=1});
c1F:=eval(a1,AG_vals[2] union {S=1.});
simplify(c1E-c1F);
c2E:=diff(c1E,eta,eta): c3E:=subs(eta=0,c2E);
c2F:=diff(c1F,eta,eta): c3F:=subs(eta=0,c2F);
c3E-c3F; # -2.*10^(-14)

You will see that the difference between c3E and c3F is definitely acceptable since Digits = 15.
If you also try this:
 

fnormal(simplify(c1E-c1F));

you will find 0. fnormal removes small terms. See the help for fnormal.

Try this:
 

plots:-implicitplot(S_TEV-S_P=0,E = 0 .. 500000, T = 0 .. 15);

@RLessard Thanks for the screenshot. That convinces me that indeed it was looking as you want it.

That I couldn't remember seeing that ever can be explained as a combination of my old age and the fact that I almost never used Greek letters as formal arguments: They are too long and what is the point?
Anyway, try this in your more recent Maple versions:
 

p:=(alpha,beta)->alpha+beta;

I'm using Maple 2021.2 on this computer and the printing of the procedure p uses Greek characters.
Because of the options operator,arrow the following is printed exactly the same:

q:=proc(alpha,beta) option operator,arrow;  alpha+beta;  end proc;

In Maple 2021.2 (at least) also these versions having a local work: as you want
 

q1:=proc(alpha,beta) option operator,arrow; local b;
   b:=7;
    b*(alpha+beta);  
end proc;
q1(1,2);
q2:=(alpha,beta)-> local b:=7; b*(alpha+beta);
q2(1,2);

In Maple 12 q2 doesn't work. It may not work in your versions because the syntax
(alpha,beta)-> local b:=7; b*(alpha+beta);
was introduced in a rather recent release.

PS. Somewhat late I paid attention to your point about saving the worksheet and opening it again.
If you save with output (I very rarely do) then indeed the Greek characters are there.

You most likely don't really like to see the steps Maple is actually taking. Take a look at this:
 

int(2/(1-x^p)^(1/p), x=0..1,method=_RETURNVERBOSE);

The answer is:
[meijerg = 2*Pi/(p*sin(Pi/p)), FAILS = (distribution, piecewise, series, o, polynomial, ln, lookup, cook, ratpoly, elliptic, elliptictrig, meijergspecial, improper, asymptotic, ftoc, ftocms, contour)]
which means that the only successful method open to int is method = meijerg.

You need 3 members in your scene. Use DEplot3d.
 

restart;
mu := 1/(365*80);
alpha := 20;
psi := 100;
varphi := 2;
K := 100;
theta := 0.5;
sigma := 0.001;

with(DEtools):
with(PDEtools):
beta := 0.3;
des := diff(S(t), t) = (psi + mu)*K - beta*S(t)*In(t) - (mu + varphi + psi)*S(t) + (theta - psi)*V(t) - psi*In(t), diff(V(t), t) = varphi*S(t) - sigma*beta*In(t)*V(t) - (mu + theta)*V(t), diff(In(t), t) = beta*S(t)*In(t) + sigma*beta*In(t)*V(t) - (mu + alpha)*In(t);

#scene = [S(t),V(t), In(t)]:
DEplot3d([des], [S(t), V(t), In(t)], t = 0 .. 100, [[S(0) = 70, V(0) = 0, In(0) = 30]], V = 0 .. 100, In = 0 .. 100, stepsize = 0.05, arrows = medium, scene = [S(t),V(t), In(t)], linecolour = sin(t*Pi/2));

The output from the first is L. I think you want NULL, because you are just making a global change on L:

f:=proc(Z)
   global L;  
   #do some processing on Z, then collect it into list L
   L:=[ op(L), Z^2];
   #print("op(L) = ",op(L)," Z=",Z); 
   NULL;  
end proc;