(F must be even.)

alpha := (1/8)*Pi:  R := 100: Wi := 3/10: H := 0: 
Eq1 := diff(F(eta), eta, eta, eta)+2*alpha*R*F(eta)*(diff(F(eta), eta))+4*alpha^2*(diff(F(eta), eta))+Wi*(8*alpha*F(eta)*(diff(F(eta), eta, eta, eta))+32*alpha^3*F(eta)*(diff(F(eta), eta)))-H*R*(cos(alpha*eta)^2*(diff(F(eta), eta))-alpha*F(eta)*sin(2*alpha*eta)) = 0:
IC1 := F(0) = 1, F(-1) = 0, F(1) = 0:
MyIC1:= F(0) = 1, D(F)(0) = 0, F(1) = 0:
MyIC2:= F(0) = 1, D(F)(0) = 0, F(-1) = 0:
smy1:=dsolve({Eq1, MyIC1}, F(eta), numeric):
smy2:=dsolve({Eq1, MyIC2}, F(eta), numeric):
p1:=plots:-odeplot(smy1):
p2:=plots:-odeplot(smy2):
plots:-display(p1,p2);

Explore could be helpful. It is not clear if a solution exists.
 

s:=dsolve({Eq1, F(0)=1, D(F)(0)=p, (D@@2)(F)(0)=q}, F(eta), numeric, parameters=[p,q]);
P:=proc(p,q)
  s(parameters=[p,q]);
  plots:-odeplot(s, -1..1)
end;
Explore( P(p,q), p=-10. .. 10., q=-10. .. 10., initialvalues=[p=0,q=0] ); 

Edit: forgot the interval in odeplot

Set
Digits:=45;
and it works. (40 is not enough.)
t0 := 0.551325384496649794326454673778951177891621325e-1

The question is why do you need Analytic for a simple equation in x
( cos(Pi*x*p)=q).

`index/defseq1` := proc(L)   local n:=op(L);  
      n^2   end proc:
a:=table(defseq1);

a[10] + a[n+1];
   100+(n+1)^2

You are not in a Boolean context, so you must use:

(t->evalb(t="."))(9);

You have used FDM in another interval.
Setting N=5 the two methods will agree.
If you keep N=10, you must set h:=1 but then fsolve cannot solve the system;  try to use DirectSearch instead.

If you want the big picture of the solutions of your system, of even if you want all the solutions in a prefered order, use Groebner bases.

For example, a system with 3 unknowns

P1(x,y,z,a)=0, P2(x,y,z,a)=0, P3(x,y,z,a)=0.

can usually be rewritten in an equivalent form:

Q1(x,a)=0, Q2(x,y,a)=0, Q3(x,y,z,a)=0.
and here the structure is clear (Pi,Qi are polynomials).
See ?Groebner
But as Carl said, we cannot treat an "abstract" case without an example, it would imply a course in Polynomial algebras, ideals etc.

Any 3D vector graphics is terrible or simply not working. The old Classic Interface is a little better but this interface is now absent for a 64 bit OS.
I have posted in the past some examples here but I cannot find them.
It seems that Maple considers this (very old!) problem as irrelevant. I use Maxima (free) for vector 3d graphics.
Here is how a simple plot
plot3d(x^2+y^2, x=-1..1,y=-1..1);
is exported as encapsulated postscript (in pdf it is exported as bitmap!):

P:=(u,v)->[u,v,-(u^2+v^2)]:
r:=(u,v,t)-> `if`(-2*u*cos(t)-sin(t)<0,0.1,1):
g:=(u,v,t)-> `if`(-2*u*cos(t)-sin(t)<0,0.1,0):
b:=g:
plots:-display(
seq( plot3d( P(u,v), u=-1..1, v=-1..1,  color=[r(u,v,t),g(u,v,t),b(u,v,t)]), t=0..2.*Pi,Pi/20.),
insequence=true);

Or, you can use Explore:

Explore( plot3d( P(u,v), u=-1..1, v=-1..1,  color=[r(u,v,t),g(u,v,t),b(u,v,t)]),  t=0..2*Pi, animate, loop,autorun) ;

If d(x,y,z) is the divergence, the integral is:

int( d(x,y,z), [z=x^2+y^2..-2*x, y=-sqrt(-2*x-x^2)..sqrt(-2*x-x^2), x=-2..0]);

Plot:
plot3d(  [x^2+y^2,  -2*x],  x=-2..0, y=-sqrt(-2*x-x^2)..sqrt(-2*x-x^2));

plot3d of
s:=[sin(phi)*cos(theta), sin(phi)*sin(theta), cos(phi)];
for phi, theta satisfying

{op(s <= ~ a), op(s >= ~ -a),  phi<=Pi, Phi>=0, theta<=2*Pi, theta>=0}

Unfortunately solve is not able to solve this system of inequalities, so you will need to do it by hand (which is not difficult; you may also use the symmetry).

The fact that
asympt(Zeta(1/2+I*y),y,6): simplify(%);
     1+2^(-1-(2*I)*y)+3^(-1/2-I*y)+2^(-1/2-I*y)

"works"  is just a mistake.  Please note the absence of the O(...) term.
Note also that
MultiSeries:-asympt(Zeta(1/2+I*y),y,6);
==> Error, (in MultiSeries:-multiseries) unable to compute multiseries

N:=1000:
f:=n ->evalf(frac(n*sqrt(2))):
g:=n ->evalf(frac(ln(n))):
X:=Vector(N,i->i):
Yf:=Vector(N,f):  Yg:=Vector(N,g):
plots:-display( Array([
plot(X,Yf, style=point), plot(X,Yg, style=point)]) );

Here is a crude version of a Frullani proc.

Fru:=proc(f,x,a,b,x0:=1)
local A,B,X,JA,JB;
B:=limit(f, x=infinity);
if type(B,`..`) or has(B,{limit,undefined}) then
  JB:=int(f,x=x0..X) assuming X>1000*(x0+1);
  B:=limit( JB/X, X=infinity)  fi;
A:=limit(f,x=0,right);
if type(A,`..`) or has(A,{limit,undefined}) then
  JA:=int(f/x^2,x=X..x0) assuming X>0,X<x0;
  A:= limit(X*JA, X=0, right) fi;
Int( (eval(f,x=a*x) - eval(f,x=b*x))/x, x=0..infinity) = (B-A)*ln(a/b)
end;

#tests

Fru( abs(sin(x)),x,2,1);
Fru( sin(sin(x))/x, x, a,b);

Edit: forgot an "x" in Int