I didn't get any empty plot when plotting the following randomly generated expression. As expected it oscillates a lot.
 

restart;
r:=rand(-350..350);
r(); # Test of r
H:=n->add((r()+r()*I)*exp(r()*I*t),i=1..n);
H(9); # Test of H
plot(abs(H(50)),t=0..2*Pi,size=[1800,default]);
plot(abs(H(500)),t=0..2*Pi,size=[1800,default]);

interface(version);

Standard Worksheet Interface, Maple 2015.2, Windows 8, December 21 2015 Build ID 1097895

(I'm using Maple 2015.2 on a 64bit Windows 10 computer).

I didn't have complete success and furthermore I found a bug.
I do have a good guess about the correct value: See at bottom.
Notice that your function P is symmetric with respect to the origin, so the two integrals A1 and A2 below ought to be equal.
As it turns out A2 comes out as a nice looking result. Regrettably it is wrong.
 

restart;
P := (x,y)->(1/2)*exp(-(1/2)*(x^2+G*y^2-2*B*x*y)/(-B^2+G))/(Pi*sqrt(-B^2+G));
P(x,y)-P(-x,-y); #Symmetry w.r.t. (0,0)
A1:=Int(Int(P(x,y),x=0..infinity),y=-infinity..0); #The first integral
## Now making a simple change of variable in A1:
PDEtools:-dchange({x=-xx,y=-yy},A1,[xx,yy]);
A2:=subs(xx=x,yy=y,%); 
#This is the same as your second integral, i.e. the integral over the second quadrant.
res1:=value(A1) assuming G>B^2; #The inner integral is computed, but not the outer.
## A numerical test:
evalf(eval(res1,{B=-1,G=2}));
evalf(eval(A1,{B=-1,G=2}));
## The simple test was passed.
## Now try also A2 (although that ought to give the same as A1)
res2:=value(A2) assuming G>B^2; #Surprisingly simple!
## The same numerical test fails:
evalf(eval(res2,{B=-1,G=2}));
evalf(eval(A2,{B=-1,G=2}));
## Converting to polar coordinates doesn't appear to help:
A1pol:=Int(Int(P(r*cos(theta),r*sin(theta))*r,r=0..infinity),theta=-Pi/2..0);
value(%) assuming G>B^2; #Worthless
evalf(eval(A1pol,{B=-1,G=2})); # Agrees with the res1 test.

It seems that the correct value is given by
 

res2p:=piecewise(B<0,res2+1/2,res2);
## Numerical evidence:
plot3d(res2p-A2,B=-1..1,G=1.1..2);

Note:
Instead of using a piecewise description of the presumed correct value we can use the two argument arctan:
 

res2p:=piecewise(B<0,res2+1/2,res2);
res3:=arctan(sqrt(-B^2+G),B)/2/Pi;
simplify(res2p-res3) assuming B>0,G>B^2;
simplify(res2p-res3) assuming B<0,G>B^2;

Here is a way. (I removed the `` around the records).
 

A:=Vector[column](4, [J_K = Record(mu = 724.901557888305, sigma = 96.7437910529146), I_W = Record(mu = 775.098442111694, sigma = 96.7437910529198), K_J = Record(mu = 785.098442111694, sigma = 96.7437910529198), D_B = Record(mu = 764.901557888305, sigma = 96.7437910529146)]);
B:=rhs~(A);
p:=sort(B,(x,y)->x:-mu>y:-mu,output=permutation);
A[p];

L:=[[TC,DB], [], [TD,JK], [IW,CM], [], [KJ,DJ]];
subs([]=NULL,L);

or use eval:
 

eval(L,[]=NULL);

I think that your limits are not covered by the numerical routines.
In the help page for int, numeric it says in the section Outline of the Numerical Integration Polyalgorithm:

For the limits of integration, the values infinity and/or -infinity are valid ...

I think that this rather indirectly must imply that your limits are not acceptable.

The value (indeed evalf(Pi) ) comes via `evalf/int` and then `evalf/int/control` from `evalf/int/improper`.
Try the direct call:

`evalf/int/improper`(1/(1+s^2),s,1-I*infinity,1+I*infinity,.5e-9,_Gquad,[]);

In line 9 of `evalf/int/improper` (under the else part after line 4 ) we are led to
2*procname(f1,x,0,infinity,1/2*eps,meth,methopts)
because  f = 1/(1+s^2) is even in s, so f1-f2 =0.
Notice that 'infinity' is hard coded into that recursive call to `evalf/int/improper`.
So basically the reason we get Pi is that 2*int(1/(1+s^2),s=0..infinity) = Pi.

Besides remembering the assignment mentioned by Kitonum, I suggest using the recommended add  instead of sum in the finite sums. Then the assignment to k will be no problem.

But I wonder why you use sum in the first place since the upper limit is the same as the lower in all cases.
This part of the code:
 

evalf(sum(S2, k = 5 .. 5)), evalf(sum(S2, k = 10 .. 10)), evalf(sum(S2, k = 20 .. 20)), evalf(sum(S2, k = 50 .. 50)), evalf(sum(S2, k = 100 .. 100)), evalf(sum(S2, k = 200 .. 200)),evalf(sum(S2, k = 500 .. 500));

could be written much easier by using seq:

evalf(seq(S2,k=[5,10,20,50,100,200,500]));

First a comment.
In order to help you better you ought to tell us also what you are actually trying to do in mathematical terms.
Is H(t) in fact supposed to be (wl/gamma1)*cos(w*t)* diff(x(t),t) ?
And is e supposed to be  int(m0*H(t),t=0..2*Pi/w)  ?
### Now just code for dsys using the parameters option in dsolve/numeric.
 

restart;
## Assigning to the parameters that won't get changed:
wl:=0.08;
alpha:=.5;
gamma0:=8.82*1e10;
mu0:=4*Pi*1e-7;   
gammA:=gamma0/(1+alpha^2);
alphaPrime:=mu0*gammA*alpha;
gamma1:=mu0*gammA;
m0:=mu0*1e5;
##
dsys:={diff(x(t), t) =alphaPrime*(z(t)^2*wl*cos(w*t)/gamma1+y(t)^2*wl*cos(w*t)/gamma1),diff(y(t),t)=-z(t)*wl*cos(w*t)-alphaPrime*y(t)*wl*cos(w*t)*x(t)/gamma1,diff(z(t),t)=y(t)*wl*cos(w*t)-alphaPrime*z(t)*wl*cos(w*t)*x(t)/gamma1,x(0)=.01,y(0)=1,z(0)=0.1};
## Using the parameters option with parameter w:
resN:=dsolve(dsys,numeric,parameters=[w],maxfun=0,range=0..1);
##
## The procedure Q sets the parameter in resN and just returns resN with the parameter set.
## Thus e.g. Q(200) would return resN with parameter w set at 200.
##
Q:=proc(w) if not w::realcons then return 'procname(_passed)' end if;
  resN(parameters=[w]); 
  resN 
end proc;
## A plot and some animations:
plots:-odeplot(Q(100),[t,x(t)],0..1,size=[800,default],refine=1);
plots:-animate(plots:-odeplot,[Q(w),[t,x(t)],0..1,refine=1],w=100..500,frames=50,size=[800,default]);
plots:-animate(plots:-odeplot,[Q(w),[t,y(t)],0..1,refine=1],w=100..500,frames=50,size=[800,default]);
plots:-animate(plots:-odeplot,[Q(w),[t,z(t)],0..1,refine=1],w=100..500,frames=50,size=[800,default]);

Here is the first animation:

Now if your answer is yes to my initial questions, then you can find e easily by including and extra unknown function Ha(t) in the system dsys as follows:

dsysH:={diff(x(t), t) =alphaPrime*(z(t)^2*wl*cos(w*t)/gamma1+y(t)^2*wl*cos(w*t)/gamma1),diff(y(t),t)=-z(t)*wl*cos(w*t)-alphaPrime*y(t)*wl*cos(w*t)*x(t)/gamma1,diff(z(t),t)=y(t)*wl*cos(w*t)-alphaPrime*z(t)*wl*cos(w*t)*x(t)/gamma1,
diff(Ha(t),t)=m0*(wl/gamma1)*cos(w*t)* diff(x(t),t),
Ha(0)=0,x(0)=.01,y(0)=1,z(0)=0.1};

Proceed as above by finding resN and Q.
To find e at some value of w, say w=200, just do Q(200)(2*Pi/200) and look at the result for Ha. I find for this example: 7.19660950992887*10^(-11).
Here is more of an automated handling using output=listprocedure in dsolve and making another procedure QH similar to Q but returning Ha(2*Pi/w):

resN:=dsolve(dsysH,numeric,parameters=[w],maxfun=0,range=0..1,output=listprocedure);
HA:=subs(resN,Ha(t)); # The procedure for computing Ha.
QH:=proc(w) if not w::realcons then return 'procname(_passed)' end if;
  resN(parameters=[w]); 
  HA(2*Pi/w) 
end proc;
## Testing QH:
QH(200);
# Plotting what I'm guessing is your e as a function of w:
plot(QH(w),w=100..10000);

If I'm right in my guesses then the relevant animation of x(t) might rather be:

plots:-animate(plots:-odeplot,[Q(w),[t,x(t)],0..2*Pi/w,refine=1],w=100..500,frames=50,size=[800,default]);

where the change is that the t-interval is now 0..2*Pi/w rather than 0..1.

There is no reason to make the problem more complicated than it already is by differentiating one of the equations as you do.
The original system sys3 can be solved for the highest derivatives in there:
 

solve(sys3, {diff(f1(x), x$4), diff(f2(x), x$4), diff(f3(x), x$6), diff(f4(x), x, x), diff(f5(x), x, x)});

So stay with sys3 except that you should make the usual change of substituting omega^2 = omega2:
 

newsys:=subs(omega^2=omega2,sys3);

Define extra_bcs as you do originally. Then do the loop:
 

for bb in extra_bcs do 
  try 
    print(bb = .1); 
    res[bb] := dsolve(newsys union {bb = .1} union bcs3, numeric, initmesh=4096,abserr=1e-2,
      approxsoln = [omega2 = 1, f1(x) = x^4*(1-x)^4, f2(x) = x^3*(1-x)^3, f4(x) = x*(1-x), f3(x) = x^5*(1-x)^5, f5(x) = x*(1-x)]) 
  catch: 
    print(lasterror) 
  end try 
end do; 

Experiment with different values of the options (initmesh, approxsoln, abserr). As they are presented here I don't get any results.
 

The second example can (should) be handled this way (where for clarity I add 7 to the action on z(x)):
 

restart;
eq := diff(y(x), x, x)+y(x) = 0;
nan1 := dsolve({eq, y(0) = 1, (D(y))(0) = -1,z(0)=0}, numeric, 
discrete_variables=[z(x)::float],events = [[y(x) = 0, [z(x) = x*y(x)+7,halt]]]);
nan1(20);

Notice that an initial condition must be given for z.

Here is a working method that selects the names of the form xxx__b:

u:=a__b+b__a+a__b^2+7+a*b+q__b;
nms:=convert(indets(u,name),list);
S:=map(convert,nms,string);
res:=StringTools:-Search~("__b",S);
S1:=[seq(`if`(res[i]<>0,S[i],NULL),i=1..nops(S))];
map(parse,S1);

You could finish with
select(has,u,%);
if that is what you want.
## Note: x__b is printed as it is because the built-in procedure print knows how to print names of that kind. I believe that no other part of the system knows any difference in principle between xxx__b and xxxyyb.

If you don't want to use the suggestion by Kitonum  g:=D[1](f); then you could do:
 

restart;
f:=(x,y)->(x^3+y^3)^(1/3);
##
g:=unapply(diff(f(x,y),x),x,y);
## Test:
g(4,5);

I find that in some situations D seems a little stupid, so often I prefer diff with unapply.
Here is an example:

restart;
f:=(x,y)->(x^3+y^3)^(1/3);
r:=(s,t)->(s*sin(t),s+56);
f(r(s,t));#OK
(f@r)(s,t);#OK
(f@r)(1,2); #OK
f(r(1,2)); #OK
D[1](f@r); #ERROR
g:=unapply(diff((f@r)(s,t),s),s,t); # Fine
g:=unapply(diff(f(r(s,t)),s),s,t); # Fine

Your worksheet works fine in my Maple 2017.2 and Maple 17.02, but works as you show in Maple 2016.2, Maple 2015.2, and Maple 18.02.
Since your f2, f3, and f4 are equations, I tried using map in Maple 2015.2:
 

f9:=map(Diff,f2,t);
f10:=map(Diff,f3,t);
f11:=map(Diff,f4,t);

That helped on the first two, but the third still had issues.
This more complicated version works on f4 as well as on f2 and f3:

evalindets(f4,function,Diff,t);

## NOTE:
I noticed a difference in the latter command if Diff is changed to %diff. 
 

evalindets(f4,function,%diff,t);

The derivative of sigma__p(t) is here shown as a dot (Newton notation) whereas it is shown with d/dt (Leibniz notation) when using Diff. This happens in Maple versions prior to 2017.
In Maple 2017.2 both (%diff and Diff) use Leibniz notation.

Here is a solution:
 

restart;
ODE:= -2*sin(1/2*theta(t))*cos(1/2*theta(t))*(diff(theta(t),t)^2-9.8000*
sin(theta(t))-(150+4*sin(1/2*theta(t))^2)*diff(theta(t),t,t))=0;
## I hope I got it right. I hate 2D math input with a passion!
##
ICS:=  theta(0) = Pi/6, D(theta)(0) = DthetaZero;
res:=dsolve({ODE,ICS},numeric,parameters=[DthetaZero],output=listprocedure);
T,T1:=op(subs(res,[theta(t),diff(theta(t),t)]));
Q:=proc(t,DTZ) option remember; res(parameters=[DTZ]); evalf([T(t)-Pi/3,T1(t)]) end proc;
q1:=proc(t,DTZ) Q(_passed)[1] end proc;
q2:=proc(t,DTZ) Q(_passed)[2] end proc;
##
L:=fsolve([q1,q2]);
res(parameters=[L[2]]);
res(L[1]); # Compares favorably with theta(t) = Pi/3 and theta'(t) = 0.
evalf(Pi/3);
## A rather flat top:
plots:-odeplot(res,[[t,theta(t)],[t,diff(theta(t),t)]],L[1]-.5..L[1]+.5,thickness=3); p1:=%:
p2:=plot([L[1],th,th=0..Pi/3],color=red,thickness=3): 
plots:-display(p1,p2);

The graphs represent theta and theta' near the t-value found.

1. Your expression f1 is (already) an equation of the form xxx = 0, so f1 = 0 makes no sense.
    (I notice that the worksheet has f1 = 0, but that the code in the question doesn't.)
2. plot doesn't plot equations, but expressions in one variable.
3. Your equation f1 contains two variables:  f and Vx. Since they are related by f1 you could try using implicitplot as in

plots:-implicitplot( f1, f = 10..1e5, Vx = 10..1e5);

Here substitute a range for Vx that you think is reasonable.
I haven't had any luck!

## Are kzp and kzm ever real?
## Yes indeed, see my reply below.

Besides returning unevaluated as in Carl's elegant version the use of Pi (as it appears in theta and theta1) may have to be handled. That could be done in the procedure at.
Here I pick rather arbitrarily an ode myself to illustrate what could be done.
I'm using a more clumsy looking version of returning unevaluated than Carl's:
 

restart;
ode:=diff(f(y),y,y)+sin(f(y))=y; #Example
nans:=dsolve({ode,f(0)=1,D(f)(0)=0},numeric);
## The important change is in the very first line:
at := proc (xx, yy) local x:=evalf(xx);
   if not type([x,yy],list(realcons)) then return 'procname(_passed)' end if;
   if 0 < x then arctan(yy/x) elif x < 0 then Pi+arctan(yy/x) elif x = 0 then (1/2)*Pi end if 
end proc;
##
theta := at(f(y)-(1/2)*Pi, diff(f(y), y));
##
plots:-odeplot(nans, [[y, theta]], y = 0 .. 6);
##
theta1:=y->at(f(y)-Pi/2,diff(f(y),y));
##
plots:-odeplot(nans, [[y, theta1(y) ]], y = 0 .. 6);