The imagenary unit should be coded in Maple as  I  (not  i ).
Here is an example of solving the equation  f(w,B,V,k)=0  for  w  for specified values  B, V, k :

gm:=V->1/sqrt(1-V^2): T:=w-k*V: S:=w*V-k:
f:=unapply(I*B*gm(V)^3*T^3+gm(V)^2*T^2-I*gm(V)^3*T*S^2-1/3*I*B*gm(V)^3*S^3*T-1/3*gm(V)^2*S^2, w,B,V,k);
F:=(B,V,k)->fsolve(f(w,B,V,k), w, complex);
F(1, 0.5, 2);

If I understand correctly, you have an equation with several parameters and you want to plot a graph that shows how the roots change depending on the change of parameters. To do this, you need to solve the equation numerically using  fsolve  command. Since the roots are not one, but several, then the number of plots should be the same as the number of roots. In the example below, we explore how the imaginary part of the roots of a quadratic equation changes depending on the change in parameters  p  and  q :

f:=(x,p,q)->x^2+p*x+q:
F:=(p,q)->fsolve(f(x,p,q), complex);
P1:=plot3d(Im('F'(p,q)[1]), p=-10..10, q=-10..10, color=grey, numpoints=10000, axes=normal):
P2:=plot3d(Im('F'(p,q)[2]), p=-10..10, q=-10..10, color=grey, numpoints=10000, axes=normal):
plots:-display(<P1 | P2>);

The graphs clearly show that the region where the imaginary part is 0, i.e. the roots are real, separated by the parabola  p^2-4*q=0  (i.e. the discriminant is 0), and the values of the complex roots are conjugate.

Edit.

Download maple17_new.mw

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