Probably this model only works adequately on a narrower time range.

For example for  t=0..10 we have:

restart;
with(plots):

s:=10:  m1:=0.02:  m2:=0.5:  m3:=4.4:  r:=0.03:  Tm:=1500:  k:=0.000024:  N:=300:    A:=1:  B:=1-(T(t)+Ti(t))/(Tm): 

dx[1] := -T*V*k*u+B*T*r-T*m1+A: 
dx[2] := T*V*k*u-Ti*m2: 
dx[3] := N*Ti*m2-V*m3: 
H := A*T-(1-u)^2+add(dx[i]*L[i], i = 1 .. 3): 
satu := -(diff(H, T)): 
dua := -(diff(H, Ti)): 
tiga := -(diff(H, V)): 
empat := diff(H, L[1]): 
lima := diff(H, L[2]): 
enam := diff(H, L[3]):

eq1 := diff(L1(t), t) = -A-(-V(t)*k*u(t)+B*r-m1)*L1(t)-u(t)*k*V(t)*L2(t): eq2 := diff(L2(t), t) = -N*m2*L3(t)+m2*L2(t): 
eq3 := diff(L3(t), t) = T(t)*k*u(t)*L1(t)-T(t)*k*u(t)*L2(t)+m3*L3(t): 
eq4 := diff(T(t), t) = -T(t)*V(t)*k*u(t)+B*T(t)*r-T(t)*m1+A: 
eq5 := diff(Ti(t), t) = T(t)*V(t)*k*u(t)-Ti(t)*m2: 
eq6 := diff(V(t), t) = N*Ti(t)*m2-V(t)*m3:

u:=t->-1/(2)*L1(t)*k*V(t)*T(t)+1/(2)*L2(t)*k*V(t)*T(t)+1: 

fcns := {L1(t), L2(t), L3(t), T(t), Ti(t), V(t)}:
 
a := dsolve({eq1, eq2, eq3, eq4, eq5, eq6, L1(10) = 0, L2(10) = 0, L3(10) = 0, T(0) = 800, Ti(0) = 0.4e-1, V(0) = 1.5}, fcns, type = numeric, method = bvp[midrich]):

odeplot(a, [[t, u(t)], [t, V(t)]], 0 .. 10, numpoints = 1000);

But for  t=0..11  we have:

restart;
with(plots):

s:=10:  m1:=0.02:  m2:=0.5:  m3:=4.4:  r:=0.03:  Tm:=1500:  k:=0.000024:  N:=300:    A:=1:  B:=1-(T(t)+Ti(t))/(Tm): 

dx[1] := -T*V*k*u+B*T*r-T*m1+A: 
dx[2] := T*V*k*u-Ti*m2: 
dx[3] := N*Ti*m2-V*m3: 
H := A*T-(1-u)^2+add(dx[i]*L[i], i = 1 .. 3): 
satu := -(diff(H, T)): 
dua := -(diff(H, Ti)): 
tiga := -(diff(H, V)): 
empat := diff(H, L[1]): 
lima := diff(H, L[2]): 
enam := diff(H, L[3]):

eq1 := diff(L1(t), t) = -A-(-V(t)*k*u(t)+B*r-m1)*L1(t)-u(t)*k*V(t)*L2(t): eq2 := diff(L2(t), t) = -N*m2*L3(t)+m2*L2(t): 
eq3 := diff(L3(t), t) = T(t)*k*u(t)*L1(t)-T(t)*k*u(t)*L2(t)+m3*L3(t): 
eq4 := diff(T(t), t) = -T(t)*V(t)*k*u(t)+B*T(t)*r-T(t)*m1+A: 
eq5 := diff(Ti(t), t) = T(t)*V(t)*k*u(t)-Ti(t)*m2: 
eq6 := diff(V(t), t) = N*Ti(t)*m2-V(t)*m3:

u:=t->-1/(2)*L1(t)*k*V(t)*T(t)+1/(2)*L2(t)*k*V(t)*T(t)+1: 

fcns := {L1(t), L2(t), L3(t), T(t), Ti(t), V(t)}:
 
a := dsolve({eq1, eq2, eq3, eq4, eq5, eq6, L1(11) = 0, L2(11) = 0, L3(11) = 0, T(0) = 800, Ti(0) = 0.4e-1, V(0) = 1.5}, fcns, type = numeric, method = bvp[midrich]):

odeplot(a, [[t, u(t)], [t, V(t)]], 0 .. 10, numpoints = 1000);

Output is:
Error, (in dsolve/numeric/bvp) initial Newton iteration is not converging
Error, (in plots/odeplot) input is not a valid dsolve/numeric solution
 

simplify((-1)^n) assuming n::even;
simplify((-1)^n) assuming n::odd;

Edit. If  nx=1 , then the result will always be 1. Therefore, Maple will not be able to show the correct answer without knowing what is  true: nx=1  or  nx=-1 . 

This can be done in many ways. Here is a simple method:

A1:= [a1, a2, a3];
A2:= [b1, b2];
[seq(seq(x/y, x=A1), y=A2)];

In the general case, instead of division, there will be a function of two variables  f(x,y).

Download collect_denominator_new.mw

This will be the case if you have an inert sums on the left. But according to your picture it is impossible to understand how the expressions on the left are obtained.

Inert sums are convenient for processing step-by-step calculations as in the code below:

Download inertsum.mw

The values of all the parameters are needed, for example:

restart;
plot3d(eval(kappa*sech(kappa*(-k^2*t+x)), [kappa=1,k=2]), x = -10 .. 10, t = -10 .. 10);

plot3d(eval((2*(-K^2+k^2))*(K*cosh(k*(-k^2*t+x))+k*cosh(K*(-K^2*t+x)))/((K+k)^2*cosh(k*(-k^2*t+x)-K*(-K^2*t+x))+(-K+k)^2*cosh(k*(-k^2*t+x)+K*(-K^2*t+x))+4*K*k), [k=2,K=3]), x = -10 .. 10, t = -10 .. 10);

plot3d(eval(2*(diff(arctan(lambda*sin(mu*((-3*lambda^2+mu^2)*t+x))*sech(lambda*((-lambda^2+3*mu^2)*t+x))/mu), x)), [lambda=0.5,mu=1]), x = -10 .. 10, t = -10 .. 10);

OP wrote  "remove 0=0 from ... if it lie in any position from the set? But  EQI__4  is a list, not a set. If we write  EQI__4  as a set and it does not contain equations in which the number on the left side is less than 0, then due to the automatic sorting of sets, the equation  0=0  will be the first and it is very easy to remove it:

 EQI__4:={omega__0/3=1/30,  -7*omega__0/6-omega__1/6=-13/60, 11*omega__0/6+5*omega__0/6+omega__2/3=47/60, 0=0, omega__1/3+5*omega__2/6=9/20};
EQI__4[2..-1];
 

The code (an example):
 

Eq:=a * q -2 * b * q + 3 * c * q +  d * q - a * w + 5 * b * w + 3.5 * c * w + 5/8 * d * w  = 9;
A:=[op(lhs(Eq))];
B:=map(T->`if`(nops([op(T)])=2,s[cat(i,ListTools:-Search(T,A))],op(1,T)*s[cat(i,ListTools:-Search(T,A))]), A);
`+`(op(B))=rhs(Eq);

The output:

Edit.

The integrand  K1  is periodic with a period of  0.2 . The variable  d1   only affects the shift of the graph along the z-axis. The integral over any interval that is a multiple of the period does not change with shifts, so this integral does not depend on  d1-variable. So I simply removed this variable from your code.

restart;
 h:=z->1-(delta2/2)*(1 + cos(2*(Pi/L1)*(z - L1))):
 K1:=((4/h(z)^4)-(sin(alpha)/Gamma2)-h(z)^2+Nb*h(z)^4):
 lambda:=Int(K1,z=0..1):
 L1:=0.2:
 alpha:=Pi/6:
plots:-display(Matrix(2,2, [seq(plot([seq(evalf(eval(lambda,Nb=j)),j=[0.1,0.2,0.3])],delta2=0.02..0.1, color=[red,blue,green], legend=[Nb=0.1,Nb=0.2,Nb=0.3], labels=[typeset(`δ1`), typeset(conjugate(`Δp`))], title=typeset("Effect of ", 'alpha', " when ", Gamma,"2=",Gamma2)), Gamma2=[10,20,30,40])]));

The  Bisection  procedure solves your problem. Formal parameters of the procedure: a  and  b  are the left and right ends of the segment on which we are looking for the root,  Eq  is an equation in the form   f(x)=0 , eps  is the maximum calculation error.

Bisection:=proc(a,b,Eq,eps)
local a1, b1, f, c;
a1:=evalf(a); b1:=evalf(b);
f:=unapply(lhs(Eq),x);
if f(a1)=0 then return a1 else if f(b1)=0 then return b1 else
if f(a1)*f(b1)>0 then error "Should be f(a)*f(b)<0" fi; fi; fi;
do
c:=(a1+b1)/2;
if f(c)=0 or c-a1<eps then return c else
if f(a1)*f(c)<0 then b1:=c else a1:=c; fi; fi;
od;
end proc:  

Example of use:

Bisection(1, 3, x^2-2=0, 10^(-6));
                                         1.414214134

restart;
implicit_func := x^3+y^3 = 9*x*y;
c := 2;
s := evalf(solve(subs(x = c, implicit_func)));
m1 := evalf(eval(implicitdiff(implicit_func, y, x), {x = c, y = s[1]})); 
m2 := evalf(eval(implicitdiff(implicit_func, y, x), {x = c, y = s[2]}));
m3 := evalf(eval(implicitdiff(implicit_func, y, x), {x = c, y = s[3]}));
Colors:=[yellow,blue,green]:
plots:-display(plot([s[1]+m1*(x-c), s[2]+m2*(x-c), s[3]+m3*(x-c)], x=-7..7, -7..7, color=Colors), plots:-implicitplot(implicit_func, x=-7..7, y=-7..7, color=red, thickness=2, gridrefine=3), seq(plottools:-disk([c,s[i]], 0.15, color=Colors[i]), i=1..3), scaling=constrained, gridlines);

In Maple 2018.2 everything works properly (I don't have maple 18). Try  eval  command instead of subs  then calculating  m1 :

restart;
imp_fun := -4*x + 10*(x^2)*(y^(-2)) + y^2 =11:
c := 2: 
s := evalf(solve(subs(x = c, imp_fun))); 
m1 := evalf(eval(implicitdiff(imp_fun, y,x), {x = c, y = s[1]}));
plots:-display(plots:-implicitplot(imp_fun, x=-5..5, y=-5..5, color=blue, gridrefine=5, scaling=constrained), plot([[c,s[1]]], style=point, symbol=solidcircle, symbolsize=17)); # The plot for a check 

In Maple there is a rule: if the argument of a function is of  float  type, then the function is completely calculated. An universal way to get a calculated expression is to use  evalf  command:

Addition. To fully solve a trigonometric equation, use solve command with allsolutions option. For all the solutions in some range add explicit  option:

restart;
solve(sin(x)=1/2, allsolutions);
solve({sin(x)=1/2, x>=0, x<2*Pi}, allsolutions, explicit);
                      (2/3*Pi)*_B1+(1/6)*Pi+2*Pi*_Z1
                        {x = (1/6)*Pi}, {x = 5*Pi*(1/6)}

In the output  _Z1  is an integer and  _B1  is equal  0  or  1

Download evalDoubtSinx_new.mw

restart;
plot([x^(1/2), x, x^2], x=0..2, 0..2, linestyle=[dot, dash, dashdot], color=red, thickness=2, scaling=constrained, legend=["1st mode","2nd mode","3nd mode"], legendstyle=[location=right], size=[400,300]);
                              

e in Maple is simply a symbol, not  2.71828... The function e^x  should be coded as  exp(x) . The error message says so:

Warning, expecting only range variable x in expression e^x-sin(x)+x^2 to be plotted but found name e

Edit. Your code will not work anyway, since the equation  f(x)=0  has no real roots. I changed something, now everything works:

We actually see that  abs(x[6]-x[5])<10^(-14)  and  f(x[6])=0.

Download Newton.mw