Since the problem is formulated in vector form the Maple solution should be given as such, I think.
Even Maple 2021 doesn't quite do that.
Here is the closest I got:
 

restart;
X:=<x1(t),x2(t)>;
Sys:=diff(X,t)=<-4,-2;6,3>.X+<2/(exp(t)+1),-3/(exp(t)+1)>; 
ics:=eval(X,t=0) = <0,ln(2)>;
 
sol1:= dsolve([ Sys, ics]); # Not on vector form
sol:=X=eval(X,sol1); # Now on vector form

The solution looks like this:

You have a classical Volterra-Lotka system.
All solutions are periodic. The orbits in phase space are closed.
Use odeplot with a numerical solution as done here.
 

restart;
sys := diff(x1(t), t) = 3*x1(t) - 2/1000*x1(t)*x2(t), diff(x2(t), t) = 6/10000*x1(t)*x2(t) - 1/2*x2(t);
ics := x1(0) = 1000, x2(0) = 500;
res:=dsolve({sys,ics},numeric);
plots:-odeplot(res,[[t,x1(t)],[t,x2(t)]],0..15,legend=[x1,x2],size=[1000,500]);
plots:-odeplot(res,[x1(t),x2(t)],0..7,frames=100);

The last is an animation of the orbit in phse space and is seen below.

I don't have Maple 13, but I have Maple 12. It gives the same error.
The error is NOT related to the way you produce the subscript, which is seen by this code:
 

restart;
H1:=v*P[r]-((1-P[r]^2*r/(P[r]^2*r-1))/r)^(1/2);
simplify(H1); ## Error in Maple 12
normal(H1);   ## OK
radnormal(H1);         ## Same error
simplify(H1,symbolic); ## Same error

The error is also present in Maple 15, 16, 17, 18, and 2015, but not in Maple 2016.

@nm For one session you can can keep the order once determined by giving LinearAlgebra:-Eigenvectors option remember or maybe option remember,system:
 

LA_EV:=subsop(3=(remember,system),eval(LinearAlgebra:-Eigenvectors)):

So you could compare the outputs of LA_EV and LinearAlgebraEigenvectors like this:
 

restart;
Sy:=1/sqrt(2)*Matrix([[0,-I,0],[I,0,-I],[0,I,0]]);
LA_EV:=subsop(3=(remember,system),eval(LinearAlgebra:-Eigenvectors)):
to 5 do LA_EV(Sy) end do;
to 5 do LinearAlgebra:-Eigenvectors(Sy) end do;

This only assures that if you happen to execute another time the previous result will be returned.
Between sessions this is no solution. You must sort.
 

Replace expand by solve or fsolve.
Here I use fsolve. Notice that all 3 roots are real (as it happens):
 

restart;
A1 := x^3 + n - 81*x - 6 = 0;
for n to 11 do
    fsolve(A1)
end do;

Since Z is taking values in the complex plane and E is therefore also, you would need a 4D coordinate system to illustrate E's dendence on Z. Maple uses color as the fourth dimension in complexplot3d. Take a look at the help page.

rho := 1/4; mu := 1/4; E := add(Z^k/GAMMA(k*rho+mu), k = 0 .. 5);
plots:-complexplot3d(E,Z=-5-5*I..5+5*I);

You simply need to give epsilon a numeric value. Then your code runs.
epsilon is present in exz, which is going to be printed.

@AHSAN In your worksheet solve gives its answer as a RootOf expression.
That expression is not at all useless.
 

restart;
pg1 := -(888888889*sigma^6*(-(2*k*sigma + 3*Q - sigma)/sigma^2)^(5/2) - 888888889*sigma^6*((k*sigma + 3*Q - 2*sigma)/sigma^2)^(5/2) + 3333333333*(-(2*k*sigma + 3*Q - sigma)/sigma^2)^(3/2)*k*sigma^5 + 3333333333*((k*sigma + 3*Q - 2*sigma)/sigma^2)^(3/2)*k*sigma^5 + 6666666666*(-(2*k*sigma + 3*Q - sigma)/sigma^2)^(3/2)*Q*sigma^4 - 3333333333*(-(2*k*sigma + 3*Q - sigma)/sigma^2)^(3/2)*sigma^5 + 6666666666*((k*sigma + 3*Q - 2*sigma)/sigma^2)^(3/2)*Q*sigma^4 - 3333333333*((k*sigma + 3*Q - 2*sigma)/sigma^2)^(3/2)*sigma^5 + 15000000000*k^2*lambda*sigma^2 + 60000000000*Q*k*lambda*sigma - 30000000000*k*lambda*sigma^2 + 60000000000*Q^2*lambda - 60000000000*Q*lambda*sigma + 15000000000*lambda*sigma^2)/(5000000000*(k*sigma + 2*Q - sigma)^2*sigma^3);
S := solve(pg1 = 0, sigma);
Sa:=allvalues(S):
numelems({Sa});
evalf(eval([Sa],{k = 2, Q = 3, lambda = 4}));
select(type,%%,realcons);

While evalc assumes that any variables present are real, convert/trig doesn't.
In fact the equation exp(I*z) = cos(z) + I*sin(z) holds for all z in the complex plane.

evalc(exp(I*x));
convert(exp(I*x),trig);

Here you will see the solution method:
 

restart;

ode:=diff(y(x),x)*tan(diff(y(x),x))+ln(cos(diff(y(x),x))) = y(x);
infolevel[dsolve]:=5:
dsolve(ode);
DEtools:-odeadvisor(ode); # Wrong

The printout is:

Methods for first order ODEs:
-> Solving 1st order ODE of high degree, 1st attempt
trying 1st order WeierstrassP solution for high degree ODE
trying 1st order WeierstrassPPrime solution for high degree ODE
trying 1st order JacobiSN solution for high degree ODE
trying 1st order ODE linearizable_by_differentiation
trying differential order: 1; missing variables
<- differential order: 1; missing "x" successful

If you try

indets(sol,name);

you get {O,O, a, b, x, y}.
If you continue with
 

indets(%,`local`);

you get  {O,O}

so the two O's are different locals.
Three O's appear in the expression, but two of those are equal.
If you try this:

oh:=op(indets(denom(op(1,rhs(sol))),`local`));
subs(oh=XX,sol);

you will see which two ones are equal.

The complex values must come from eq3, where you have (1 + delta*theta(eta))^n and where later n = 0.2.
Just try

(-2)^0.2; # answer: 0.9293164906 + 0.6751879524*I (the principal 5th root.).
You could then think of using surd (see the help). That turns out not to work in dsolve/numeric/bvp because the round disappears and you get surd( ... , 5. ), i.e a float instead of an integer as the second argument.
You can use eq3_new instead of eq3:

eq3_new:=subs((1 + delta*theta(eta))^n=abs(1 + delta*theta(eta))^n*signum(1 + delta*theta(eta)),eq3);

The problem is illustrated here:

surd(-2, round(1/0.2)); #Here it works because the integer 5 remains an integer
evalf(%);
### The alternative
abs(-2)^0.2*signum(-2);

You run into singularity problems though, so there are obviously other problems with your system.
With these extra arguments to dsolve
method=bvp[midrich],initmesh=1024
you get the error:

Error, (in dsolve/numeric/bvp) matrix is singular

The exp term in eq3 gives problems here:
 

exp(-E/(1 + delta*theta(eta)));
####
plot(exp(-10/(1 + 10*theta)),theta=-0.5..0.5,thickness=3);

This is really just an observation, but it appears that the mere presence of sqrt(-x) or sqrt(-y) makes the difference:

restart;
A:=(x+y-2*sqrt(x*y))+sqrt(-x);
simplify(A) assuming negative;
restart;
A:=(x+y-2*sqrt(x*y))+sqrt(-x);
simplify(A) assuming positive;

Thus this trick:

restart;
# Example 2

B:=x+y-2*sqrt(x*y);

simplify([B]+sqrt(-x)) assuming negative;

I notice that solving for diff(y(x),x) and applying dsolve works (fast):
 

ode1,ode2:=solve(x^3*diff(y(x),x)^2+x*diff(y(x),x)-y(x) = 0,{diff(y(x),x)});
dsolve(ode1);
dsolve(ode2);

restart;
ode1:=diff(u(x), x, x) = u(x)*sin(x);
ode2:=diff(ode1,x);
bcs:=u(0)=2,u(1)=1,(D@@2)(u)(1)=sin(1);
Digits:=10: # The default
res:=dsolve({ode2,bcs},numeric,abserr=1e-13);
plots:-odeplot(res,[x,diff(u(x),x,x)-u(x)*sin(x)],0..1,caption=typeset(diff(u(x),x,x)-u(x)*sin(x)));

A comment or rather a rhetorical question:

If you were to design an algorithm for dsolve/numeric (bvp or ivp) that returned the highest derivative too, what would you make it return other than the right hand side?
#####
Looking at what `dsolve/numeric` does to an ivp like this:
 

restart;
ode1:=diff(u(x), x, x) = u(x)*sin(x);
ode2:=v(x)=diff(u(x),x,x);
ics:=u(0)=2,D(u)(0)=0;
#debug(`dsolve/numeric`);
res:=dsolve({ode1,ode2,ics},numeric,abserr=1e-13,relerr=1e-10);
plots:-odeplot(res,[x,v(x)-u(x)*sin(x)],0..1,caption=typeset(v(x)-u(x)*sin(x)));

If you remove the # then the debugging output reveals what dsolve/numeric does here:
 

INFO["proc"] := proc(N, X, Y, YP) local Z1; Z1 := Y[1]*sin(X); Y[3] := Z1; if N < 1 then return 0; end if; YP[2] := Z1; YP[1] := Y[2]; 0; end proc;

As you see the right hand side is kept in the local Z1 and Y[3] (i.e. v) gets that value as does YP[2] i.e. u''.
This simply means that the right and side is returned for v.