You may have expected something else than what comes out of the following slightly corrected version, if so let us know:
 

restart;
n:=5;
R:=Matrix(5,n);

for i from 1 to n do 
  R[..,i]:=sin(i*x);
end do;
R;
for i from 1 to 5 do
  #i:=i; # This line doesn't do anything useful, but isn't harmful
  x:=i;
  R[i,..] := sin(x);
end do;
R;

dsolve assumes that c and g are constants, so the first 4 solutions are constant solutions and 2 of these are imaginary when c and g are positive, which you assume.

I assume that you aren't interested in constant solutions. Thus the fifth solution is the one you want.

I don't think you can convert the JacobiSN function into a trigonometric function.

Here is an attempt to explain the situation when using assume.MaplePrimes22-09-29_assume.mw

Download MaplePrimes22-09-29_assume.mw

@acer I notice in Carl's second example that unwith followed by forget will restore the default diff.
That doesn't explain what's going on, but at least a restart is not necessary.

restart;
expr:=A/B+C;
#combine(expr);
normal(expr);

In Maple gamma is Euler's constant (see ?gamma), which is approximately equal to 0.5772156649. Try evalf(gamma).
Either of the remedies proposed by vv or tomleslie will work, however, in your case since either one is unaffected by that fact as only alpha*gamma appears. I also pointed that out in my answer to your previous question dealing with the same system. (https://mapleprimes.com/questions/234300-How-Can-I-Fix-Error-in-Dsolvenumericprocessinput ).
In general if you want to use gamma as just another name like alpha and beta, you can declare gamma local at the top of your worksheet as in:

restart;
local gamma;
## The rest ......

This only works in recent Maple versions. The one you use, Maple 2018, is included.
 

I deleted all the output that should not have been copied as Tom Leslie says.

Besides I corrected some typos and errors. So in that process I may have made my own errors.
The following appears to be what you want. Notice in particular that alpha*gamma never appears in your original equations, but alpha*m does. I changed that product to a single name: am.
 

restart; 

with(plots):  

U := [W(t), X(t), Y(t), Z(t)]; 
 f := (W, X, Y, Z) -> A-mu*W-delta1*Y-delta2*Z; 
 g := (W, X, Y, Z) -> beta*(W-X-Y-Z)*X-mu*X-X; 

h := (W, X, Y, Z) -> am*X-(mu+delta1)*Y; 

i := (W, X, Y, Z) ->(-am+1)*X-(mu+delta2)*Z; 

sys1 := diff(U, t) =~ [f(op(U)), g(op(U)), h(op(U)), i(op(U))]; 

ic1 := W(0) = W__0, X(0) = X__0, Y(0) = Y__0, Z(0) = Z__0;

equil := solve({f(w0, x0, y0, z0) = 0, g(w0, x0, y0, z0) = 0, h(w0, x0, y0, z0) = 0, i(w0, x0, y0, z0) = 0}, {w0, x0, y0, z0}); 

param1 := {A = 1, beta = 0.9e-3, delta1 = 0.139e-2, delta2 = 0.134e-2, mu = 0.132e-2, W__0 = 750, X__0 = .1, Y__0 = .2, Z__0 = .6, am = 0.4e-2}; 

param2 := {A = 1, beta = 0.9e-3, delta1 = 0.139e-2, delta2 = 0.134e-2, mu = 0.132e-2, W__0 = 755, X__0 = .2, Y__0 = .4, Z__0 = .9, am = 0.4e-2}; 

param3 := {A = 1, beta = 0.9e-3, delta1 = 0.139e-2, delta2 = 0.134e-2, mu = 0.132e-2, W__0 = 760, X__0 = .3, Y__0 = .6, Z__0 = 1, am = 0.4e-2}; 

param4 := {A = 1, beta = 0.9e-3, delta1 = 0.139e-2, delta2 = 0.134e-2, mu = 0.132e-2, W__0 = 765, X__0 = .4, Y__0 = .7, Z__0 = 1.4, am = 0.4e-2}; 

sol1 := dsolve(eval({ic1, op(sys1)}, param1), U, numeric); 

odeplot(sol1,[t,X(t)],0..30);
for j to 4 do
   p[j]:=odeplot(sol1,[t,U[j]],0..30);
end do:
display(Array([seq(p[j],j=1..4)]));

The plot of X(t):
MaplePrimes22-06-06_odesys.mw

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]