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);

Here is something to get you going:
 

restart;
Digits:=15:
## Assigning only to the parameters you don't want to change:
omega:=2*Pi*11.2e06: tau:=15: alpha:=1: gamma1:=0.06:
##  
sys:={ diff(x(t),t)=v(t)/omega, diff(v(t),t)=G*epsilon*(2*BesselJ(1, v(t-tau))+( alpha*v(t)-(3/4)*gamma1*v(t)^3))-omega*x(t)-v(t)*epsilon};
ics:={x(0)=0,v(0)=1}; # Rather arbitrary
## In dsolve/delay this will imply that x(t) = 0 and v(t) = 1 for t < 0.
params:=[G=2.5,epsilon=0.014]; #Example
## Solving just that example:
res:=dsolve(eval(sys,params) union ics,numeric);
plots:-odeplot(res,[10^8*x(t),v(t)],100..107);
plots:-odeplot(res,[10^8*x(t),v(t)],0..100,numpoints=10000,frames=25);
### Now making a procedure resf that accepts G,epsilon values.
### If not fed with actual numbers it returns unevaluated.
### That can be very convenient.
resf:=proc(G1,eps) if not type([G1,eps],list(realcons)) then return 'procname(_passed)' end if;
      dsolve(eval(sys,{G=G1,epsilon=eps}) union ics,numeric)
end proc;
## This call returns unevaluted since G is just a name:
resf(G, 0.014);
## This call returns the result res found directly above:
resf(2.5,0.014);
## Here is an animation example, where animation is in G:
plots:-animate(plots:-odeplot,[resf(G,0.014),[x(t),v(t)],400..407],G=1..7);

I would have used the parameters option in dsolve/numeric/ivp if it had been available. I forgot that it isn't.  I suppose it is documented somewhere. Therefore the need for the procedure resf.
## Actually, I don't believe that it is documented. I once filed an SCR about incomprehensible error messages when trying to use parameters as in this example:

restart;
dde:=diff(y(t),t)=y(t-1)+sin(omega*t);
res:=dsolve({dde,y(0)=1},numeric,parameters=[omega]);
res(parameters=[omega=2]);
res(3); # ERROR, message is
## 
##Error, (in res) dependent variable index out of range
##
## Try the same one more time
res(3); # ERROR, message is
##
## Error, (in res) argument #1, the current solution, delay variable specification disagrees with problem

### NOTE:
### If for your purposes you need to make one or both of the initial values x(0) and y(0) parameters, then you can use a modified version of resf as in this example where I make v(0) a parameter.
 

ics2:={x(0)=0,v(0)=v0};
resf2:=proc(G1,eps,v00) if not type([G1,eps,v00],list(realcons)) then return 'procname(_passed)' end if;
      dsolve(eval(sys union ics2,{G=G1,epsilon=eps,v0=v00}),numeric)
end proc;
resf2(G,0.014,1); #test as for resf
## Example runs using v0 = .1 and v0=5, respectively:
plots:-odeplot(resf2(2.4,0.014,.1),[x(t),v(t)],0..7,numpoints=10000,frames=10);
plots:-odeplot(resf2(2.4,0.014,5),[x(t),v(t)],0..7,numpoints=10000,frames=10);

From the last two animations you see that the orbit is on its way out in the first one, but in the second one it is on its way in.

The two boundary conditions IC1[1] and IC1[3] are equivalent, i.e. they say the same thing.

There is another problem: Confusion w.r.t. the use of the name a.
You begin by assigning to a[1], a[2], ..., a[6]. This means that now a is a table (implicitly defined by your assignments).
Then you also assign to a itself. You have a:=1. Thus the above assignments are forgotten.
Then you do
HA := [a[1], a[2], a[3]];
which now will mean 1[1], 1[2], 1[3], which clearly makes no sense. To see that try:
HA;
So don't ever use a name (here a) and its indexed versions (as in a[1] ) in the same session (unless of course by a you do mean that table I mentioned).
 

This is not an attempt to derive your solution from the odes, but here I just answer your question:
"How to create a graphic output of movement according to equation x(t) similar to graph mentioned below?"

So here is the code for producing graphical representations of the solution you give.
 

restart;
sol:=(A-n*C)*cos(omega*t)+(-1)^n*C/2;
params:={m=1.0, k=0.1, g=10.0, mu=0.003}: #Just an example
C:=2*mu*m*g/k;
omega:=sqrt(k/m);
n:=floor(omega*t/Pi);
plot(eval(sol,params union {A=1}),t=0..250, size=[800,default],caption="Initial amplitude A = 1" );
plot(eval(sol,params union {A=10}),t=0..250, size=[800,default],caption="Initial amplitude A = 10" );
## This is an animation with A as animation parameter:
plots:-animate(plot,[eval(sol,params),t=0..250],A=1..10,size=[800,default]);

The animation:

Once you have your symbolic work done you can compare it with the numerical solution using events:
 

restart;
params:={m=1,k=.1,mu=.0003,g=10};
ode:=m*diff(x(t),t,t)=-k*x(t)-mu*m*g*(2*b(t)-1);
resE:=dsolve({eval(ode,params),x(0)=1,D(x)(0)=0,b(0)=0},numeric,discrete_variables=[b(t)::boolean],events=[[diff(x(t),t)=0,b(t)=1-b(t)]]);
plots:-odeplot(resE,[t,x(t)],0..100);

The event [diff(x(t),t)=0,b(t)=1-b(t)] means that when it happens that diff(x(t),t) = 0 then b(t) is changed from its previous value (0 or 1) to the opposite, i.e. 1 or 0. That makes the factor (2*b(t)-1) change from -1 or 1 to the opposite 1 or -1.
The graph is:

If you try the range 0..200 you will discover that the events version (as I gave it) misrepresents the solution. Compare with this version which doesn't use events:
 

restart;
ode:=m*diff(x(t),t,t)=-k*x(t)-mu*m*g*signum(diff(x(t),t));
params:={m=1,k=.1,mu=.0003,g=10};
res:=dsolve({eval(ode,params),x(0)=1,D(x)(0)=0},numeric,maxfun=0);
plots:-odeplot(res,[t,x(t)],0..200);

The graph on the range 0..200:

You will notice that the solution appears to become constant at about t = 170. If this is what really happens then the event as given in the first version needs modification.