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.

Clearly a bug. And an old one at that. It appears in Maple 8 and Maple 12 as well as in Maple 2017.2.

The problem may be located in line 32 of `singular/singular` which reads:
assume(new_Z,'posint');

## I shall submit an SCR (aka a bug report).

Take a look at the help page

?Initially Known Names

You will see that gamma(0), gamma(1), ..., gamma(5), ... are known as Stieltjes constants. Furthermore gamma = gamma(0) is Euler's constant.
That is the problem, I'm sure. The unprotecting you are doing doesn't change that unless you follow it up by assigning gamma to some other name e..g. gamma:=g;
But of course that is silly: just use g to begin with instead of gamma (or h as Tom Leslie is using).

An alternative in recent versions of Maple is to place right after restart the local declaration:
local gamma;
## Omit the unprotect gamma part: it serves no purpose.
#### To see the numerical values of gamma(n) for n=0,1,2... you can try
seq(evalf(gamma(n)),n=0..5);
## To see the code for this, try
showstat(`evalf/gamma`);

The mathematical constant in question is in Maple spelled Pi, not pi. The latter is just another symbol with no meaning, but it prints the same as Pi does.
cos(Pi/2)^2;

will obviously just evaluate to 0.

You can use a try statement as in
 

Q1 := []; 
for jj to 100 do 
  simplfloat := rand(-1.0 .. 1.0); 
  h := simplfloat(); 
  a := abs(amin+(amax-amin)*h); 
  b := abs(bmin+(bmax-bmin)*h); 
  try 
    solut := dsolve([eq, cis], numeric, maxfun = 0, output = procedurelist); 
    sd := i-> abs(rhs(solut(i)[2])*rhs(solut(i)[2])); 
    eng := Simpson(sd, 5, 6, 10) 
  catch "cannot evaluate the solution": 
    printf(StringTools:-FormatMessage(lastexception[2 .. -1]));
    printf("\n"); 
    eng := FAIL 
  end try; 
  if eng = FAIL then next else Q1 := [op(Q1), [jj, a, b, h, eng]] end if 
end do:

But in my attempts with the loop going to 100 I never get any successes. So the question is: Must a singularity occur before the end of the interval?
Since what seems to happen is that x(t) decreases to zero thus creating a problem with 1/x(t), we may rephrase the question to: Must x(t) decrease to zero under the given conditions?
### Actually, I think the answer is yes, it must.
Rewriting the ode as a first order system using the usual diff(x(t),t) = y(t), we get:
 

sys:=[diff(x(t), t) = y(t), diff(y(t), t) = -(b*y(t)+sin(x(t))*y(t)+2/x(t))/a];

A small amount of phase plane analysis making use of the vertical isocline y = 0 and the horizontal isocline y = -2/x/(b+sin(x)) and the fact that a and b are positive with b > 1 will show that if the solution starts at (x,y) = (0.25, 0) (i.e. on the vertical isocline) then it will not be able to cross the horizontal isocline, which is situated in the 4th quadrant. Thus the only possibility is that x(t) -> 0 and y(t) -> -infinity as t increases to the right end of the (finite) interval of definition of the solution.
Illustration here for a = .1828083187 and b = 22.45274868.
The red line is the vertical isocline, the blue the horizontal isocline. The trajectory must stay between these thoughout its life (for t >=0 ).

Using DEplot:

Produced by:
 

plots:-display(
   DEtools:-DEplot(eval(sys,{a = .1828083187, b = 22.45274868}),[x(t),y(t)],t=0..0.38,[[x(0)=0.25,y(0)=0]],linecolor=maroon,color=gray),
   plot(-2/x/(22.45274868+sin(x)),x=0..0.25,color=blue),
   plot(0,x=0..0.25,color=red,thickness=3),
   view=-2.5..0
);

Let us consider the simplified example sol:=1/x/sin(x) to understand what is going on.
If you look at the code for algsubs:
showstat(algsubs);
you will find that you are going to `algsubs/algsubs`  with the call

`algsubs/algsubs`(1/x/sin(x),x,t,[x],remainder);

Now look at the code for `algsubs/algsubs`:
showstat(`algsubs/algsubs`); #Very short
You will see that you are going to `algsubs/recurse` with the call

`algsubs/recurse`(1/x/sin(x),x,t,[x],remainder);

Now looking at the code for `algsubs/recurse`:
showstat(`algsubs/recurse`);
you will that since type(1/x/sin(x),`*`) = true, you will be in line 4 with
map(`algsubs/recurse`, 1/x/sin(x), x, t,[x], remainder);
Thus you will handle both of these
`algsubs/recurse`, 1/x, x, t,[x], remainder);
`algsubs/recurse`, 1/sin(x), x, t,[x], remainder);
The first one will be handled by `algsubs/recurse` in line 3 since
type(1/x,'ratpoly'(anything,[x])) = true
so 1/x is returned for that factor.
As for 1/sin(x) that will be handled in line 7 of `algsubs/recurse` with the call

map(`algsubs/algsubs`,1/sin(x),x,t,[x],remainder);
resulting in 1/sin(t).
Thus altogether 1/x/sin(t).
So the crucial line is line 3 in `algsubs/recurse`:
elif type(f,('ratpoly')(anything,vars)) then f

As Kitonum says, this is not a bug. In my view because it is so blatant: It must have been intended.
As Kitonum also points out, there is this note in the help page:

"Note: The algsubs command currently works only with integer exponents."

I find that rather obscure: Integer exponents of what? And isn't -1 an integer?

## My own note: I never use algsubs, only subs and eval.

In the following I assume that the syntax errors (missing multiplication signs) have been corrected.
eq2 and eq3 are linear homogeneous odes and are subject to homogeneous boundary conditions.
So clearly psi = w = 0 (identically) is a solution. Is it the only one?
## I have looked at that question now. It is the only solution.
I have placed the proof at the end of this.
##
But let us look at that one.
So we have:
 

restart;
eq1 := (diff(psi(x), x))^2+(diff(u(x), x)+(8*(1/2))*(diff(w(x), x))^2)*(diff(psi(x), x))^2+3*(diff(w(x), x, x))+5*(diff(w(x), x, x))*(diff(psi(x), x))-7*(diff(u(x), x, x, x)+(8*(1/2))*(diff(w(x), x, x))^2+(3/2)*(diff(w(x), x, x, x))*(diff(w(x), x)))+3 = p;
eq2 := (51-31)*(diff(psi(x), x, x))+(52-2)*(diff(w(x), x, x, x))+8*(diff(psi(x), x, x, x, x))-7*(diff(w(x), x)-psi(x)) = 0;
eq3 := 4*(diff(w(x), x, x)-(diff(psi(x), x)))+(23+11)*(diff(psi(x), x, x, x))+(14+12)*(diff(w(x), x, x, x, x)) = 0;
bcs:={psi(0) = 0, psi(1) = 0, u(0) = 0, u(1) = 0, w(0) = 0, w(1) = 0, D(u)(1) = 0, D(psi)(0) = 0, D(psi)(1) = 0, D(w)(0) = 0, D(w)(1) = 0};
bcs1,bcs23:=selectremove(has,bcs,u);
eq10:=eval(eq1,{psi=0,w=0});
sol0:=dsolve({eq10} union bcs1);
plots:-animate(plot,[rhs(sol0),x=0..1],p=-50..50);

So we see that for each p there is a solution for u, psi, w where psi = w = 0, and u<>0 unless p=3:
 

solve(identity(rhs(sol0),x),p);
eval(sol0,p=3);

So you may ask, what happens numerically? Well, begin by plotting w and psi in my previous answer, where I only plotted u.
You will find that psi = w = 0.

### I place here the proof that eq2, eq3 with bcs23 only has the trivial solution psi = 0, w = 0.
 

restart;
eq1 := (diff(psi(x), x))^2+(diff(u(x), x)+(8*(1/2))*(diff(w(x), x))^2)*(diff(psi(x), x))^2+3*(diff(w(x), x, x))+5*(diff(w(x), x, x))*(diff(psi(x), x))-7*(diff(u(x), x, x, x)+(8*(1/2))*(diff(w(x), x, x))^2+(3/2)*(diff(w(x), x, x, x))*(diff(w(x), x)))+3 = p;
eq2 := (51-31)*(diff(psi(x), x, x))+(52-2)*(diff(w(x), x, x, x))+8*(diff(psi(x), x, x, x, x))-7*(diff(w(x), x)-psi(x)) = 0;
eq3 := 4*(diff(w(x), x, x)-(diff(psi(x), x)))+(23+11)*(diff(psi(x), x, x, x))+(14+12)*(diff(w(x), x, x, x, x)) = 0;
bcs:={psi(0) = 0, psi(1) = 0, u(0) = 0, u(1) = 0, w(0) = 0, w(1) = 0, D(u)(1) = 0, D(psi)(0) = 0, D(psi)(1) = 0, D(w)(0) = 0, D(w)(1) = 0};
bcs1,bcs23:=selectremove(has,bcs,u);
### Now considering only eq2 and eq3 with bcs23
sys23:={eq2,eq3};
bcs23;
## Splitting in conditions at x = 0 and x = 1:
bcs23_0,bcs23_1:=selectremove(hastype,bcs23,anyfunc(0));
## Defining initial conditions at x = 0 using names for unknown values:
ics23:=bcs23_0 union {(D@@2)(w)(0)=w2,(D@@3)(w)(0)=w3, (D@@2)(psi)(0)=p2,(D@@3)(psi)(0)=p3};
nops(ics23); # 8 as there should be
## We now have a standard initial value problem at 0 for sys23 with ics23.
## This has a unique solution for any w2, w3, p2, p3:
sol:=dsolve(sys23 union ics23): # This fills a lot, therefore the colon.
evalf(sol);# Just to see a smaller output, but will NOT be used
## Now imposing the conditions at x = 1 on the solution:
EQ1:=eval[recurse](sol,{x=1} union bcs23_1): 
evalf(EQ1); # Just having a look at a smaller version.Will NOT be used. 
## Some of the conditions at x=1 involve the first derivatives:
EQ2:=eval[recurse](convert(diff(sol,x),D),{x=1} union bcs23_1):
evalf(EQ2);# Just having a look at a smaller version.Will NOT be used.
indets({EQ1,EQ2},name);
## Now rewriting the linear system EQ1, EQ2 as a matrix equation: 
M,zz:=LinearAlgebra:-GenerateMatrix(EQ1 union EQ2,[p2,p3,w2,w3]):
zz;
evalf(M); # Just having a look at a smaller version.Will NOT be used.
## The determinant of M:
DT:=LinearAlgebra:-Determinant(M):
evalf(DT); # Apparently different from zero
evalf[50](DT); # Confirmed
is(DT>0); # Confirmed
## Thus the system M.< p2,p3,w2,w3 > = <0,0,0,0> has a unique solution.
## Since < p2,p3,w2,w3 > = <0,0,0,0> is a solution, it is the one.

My guess is that it is a simple bug. Try this:

with(ScientificConstants);
debug(Element);
GetValue(Element('H',density(gas)));
GetValue(Element('Si',density));
showstat(Element);

You will see that in line 47 of the procedure Element:
showstat(Element, 47);
ScientificConstants:-numTabElem[an]:-props[yp]:-value(`if`(withPs,params,NULL))
in the case of Hydrogen produces the number 0.088 because withPs = true, whereas for Silicium NULL is produced because withPs=false.
So my guess is that for gases the numbers produced at this stage are actually meant to be in kg/m^3, but are incorrectly treated as if in g/cm^3 and thus converted to kg/m^3 later (by GetValue in its line 42).
A closer look at line 42 in the procedure GetValue shows that in
convert(valueX,'system',ScientificConstants:-numTabElem[op(x)[1]]:-props[op(0,op(x)[2])]:-units,sys)
this part:
ScientificConstants:-numTabElem[op(x)[1]]:-props[op(0,op(x)[2])]:-units
evaluates to kg/L and sys evaluates to SI. valueX has the value 0.088.
Now kg/L is these days exactly the same as g/cm^3.
The numbers for gases must have been entered in SI units and are here are mistakenly treated as being in kg/L.

I shall submit an SCR with a link to this page.
PS. In preparing for the SCR I discovered that the density in SI units for Hydrogen in Maple 12 and Maple 15 is given as 0.000088. This changed to 88.0 in Maple 16 and all later versions.

The error message tells the story: You have a system of total order 4+4+3 = 11 and you have an unassigned parameter p.
If your intention is to determine p at the same time you solve, you must have one more boundary condition. That makes it a total of 12.

You are using the indexed variables a[i] and b[i], but also the unrelated (your intention, I assume) a and b.
That is not a good idea. As soon as you do e.g.
a[1] := 0;
a table by the name 'a' is implicitly created.
Thus it is potentially catastrophic to define as you do
a[-1] := -12*mu/(a+b);
Also a and b appear in your pde.
The solution is simple: Either change the names a and b to e.g. aa and bb (or whatever) in the pde and in
a[-1] := -12*mu/(aa+bb);
or don't use indexed variables. Use other names for a[0], a[-1], etc.
### To see the catastrophy in a simple session try this:
 

restart;
a[-1] := -12*mu/(a+b);
eval(a); #table with name a and index -1 refers to the whole table itself
simplify(a); # Not surprising!

Yes, this is a bug. You are doing nothing wrong.
The two cases are treated by separation of variables: Try
 

pdsolve(PDE1);
pdsolve(PDE2);

Both of these results are correct, just different: _c[1] in the first corresponds to c[1] - a in the second. _c[1] is an arbitrary constant at this stage.
You could go further and ask for solutions of the odes:

pdsolve(PDE1,build);
pdsolve(PDE2,build);

The solutions are on different forms for sure, but they are still solutions of their respective pdes.
So whatever goes wrong happens later.
You can if you wish try your commands after the command
debug(pdsolve);
##
I shall submit a bug report with a reference to this question.

Your problem is that this command (solving for the highest derivatives):
 

solve(newsys2,{diff(g1(y), y$5),diff(g2(y), y$4),diff(g3(y), y$6),diff(g4(y), y, y),diff(g5(y), y, y)});

has no solution.
One of the equations (in the numbering on my machine it is newsys2[4]) contains none of the highest derivatives, i.e. none of
diff(g1(y), y$5),diff(g2(y), y$4),diff(g3(y), y$6),diff(g4(y), y, y),diff(g5(y), y, y).
There is a simple solution: Differentiate that equation w.r.t. y and use that instead of newsys2[4].
Thus your new system with the highest derivatives solved for is:

SYS:=solve({newsys2[1],newsys2[2],newsys2[3], diff(newsys2[4], y), newsys2[5]}, {diff(g1(y), y$5),diff(g2(y), y$4),diff(g3(y), y$6),diff(g4(y), y, y),diff(g5(y), y, y)});

Remember though that since diff(newsys2[4],y)=0 we have that  newsys2[4] = constant, thus demanding that newsys2[4] be satisfied at one of the boundary points must be added. Thus either add this boundary condition:
 

eval[recurse](convert(newsys2[4],D),{y=0} union bcs3);

or this one
 

eval[recurse](convert(newsys2[4],D),{y=1} union bcs3);

## A comment on the boundary conditions.
You have a system that when converted to a first order system has 19 equations. You have an undetermined constant omega3, so you need 19+1 = 20 boundary conditions.
bcs3 has 18 conditions and one of the two mentioned above (from newsys2[4] at 0 or 1) makes a total of 19. You need one more condition. The extra condition can involve the following quantities:
 

N:=[5,4,6,2,2];
extra_bcs := {seq(seq(seq((D@@i)(cat(g, j))(a), i = 0 .. N[j]-1),j = 1 .. nops(N)), a = 0 .. 1)};

Since neither the conditions in bcs3 nor in the newly added one are simple i.e. of the form (D@@i)(gj)(a) = constant (a=0 or 1) it doesn't help much to define extra_bcs with an added minus bcs3.