Set initmesh higher than default, which is just 8.
You don't have to increase it much, 16 will do:
numsol1 := dsolve({BCSforNum1, BCSforNum2, ODEforNum1, ODEforNum2}, numeric, output = listprocedure, initmesh = 16);
plots:-odeplot(numsol1,[eta,u(eta)]);

You can use nops:

L:=[a+b+c,a+b,a+b+c+d,a];
nops~(L); # Notice the tilde
map(nops,L); # Alternative

parse("123");

A warning that output exceeds expression length just means that the result is not displayed on the screen. It does NOT mean that the result is not computed and cannot be used.
But as pointed out by Christopher you can change the SHOWN length of output in Tools/ Options/Precision.

You had a couple of errors. I corrected two lines:
 

dsys1[i1, i2] := eval({Eq1[i1, i2], Eq2[i1, i2], IC[i1, i2]}, m1 = 5); 
dsol1[i1, i2] := dsolve(dsys1[i1, i2], numeric, output = listprocedure, range = 0 .. L, maxmesh = 512);

I corrected a couple of errors of syntax. Be aware in particular that in Maple case matters: thus 'domain' is not the same as 'DOMAIN'.
 

restart; 
#with(plots): #Not needed
with(DEtools): 
##
ode1 := diff(x(t), t) = y(t);
ode2 := diff(y(t), t) = -x(t)+(1-x(t)^2)*y(t);
MODEL := {ode1, ode2}: 
VARS := {x(t), y(t)}; 
domain := t = 0 .. 20; #You are using lower case letters here (OK, but stay with them)
RANGE := x = -3 .. 3, y = -3 .. 3; 
COLORS := [yellow, green];
IC1 := [x(0) = .5, y(0) = .25]; 
IC2 := [x(0) = 2.5, y(0) = 3]; 

DEplot(MODEL, VARS, domain, RANGE, [IC1, IC2], stepsize = .1, arrows = THIN, linecolor = COLORS); 

You cannot have boundary conditions involving derivatives of the unknowns to more than their respective orders minus 1.
Your system has the orders: [[f1, 6], [f2, 6], [f3, 8], [f4, 4], [f5, 4]];
Thus there is a problem with bcs, which has:
select(has,bcs,(D@@8)(f3)); ## NOT empty: should be!
select(has,bcs,(D@@4)(f4));  ## same problem

Your code works for me in Maple 2016. I executed your module (received messages of two variables being implicitly declared local). Then I  did:
 

with(MonPackageFonctions);
GammaNum(.45);
psiNum(.45);

which gave as output 3.141592654-2.115788962*I and -1.369468047*I, respectively.
I assume that your r, l, t, x are globals too?
You don't have to declare GeometricData global in the procedures. Doing it in the module should suffice.
## I checked in a considerably older version, Maple 15: Same result.

I assume that your example is a toy example for a more interesting and complicated situation.
So I will show you how to use dsolve/numeric with a 'known' user defined function.
Since your example is so simple it can also be solved without the numeric option.

myfun := proc (x) local output;
if not x::numeric then return 'procname(x)' end if;
output := 4*x^2;
output
end proc;
de := diff(y(x), x)+myfun(x)*y(x) = 0.;
sol:=dsolve({de,y(0)=1}); #A formula in terms of myfun
##Now numerically
res:=dsolve({de,y(0)=1},numeric,known=myfun);
plots:-odeplot(res,[x,y(x)],0..2);
## sol can also be plotted. The integral is computed numerically because the procedure
## myfun(x) returns unevaluated if x is not of type numeric.
plot(rhs(sol),x=0..2);

It might be illuminating to see the values of x that are accessed in myfun.
To see that you can add a print statement to myfun.
 

myfun := proc (x) local output;
if not x::numeric then return 'procname(x)' end if;
userinfo(2,procname,printf("x = %f\n",x));
output := 4*x^2;
output
end proc;

Since I have wrapped the print statement in userinfo, printing is only done if infolevel[myfun] is set to at least 2 (the first argument in userinfo).
So put the line
infolevel[myfun]:=2:
before e.g. the line plot(rhs(sol),x=0..2) to see that numerical integration is indeed taking place.

Have a look at these two versions both using the sine curve.
The first represents the curve in a parametrized form, but takes sin(x) first instead of x.
The second takes an already made plot and interchanges x and y:
 

plot([sin(x),x,x=-Pi..Pi],labels=[y,x]); #A parametric version
###
p:=plot(sin(x),x=-Pi..Pi);
tr:=plottools:-transform((x,y)->[y,x]):
tr(p);
plots:-display(tr(p),labels=[y,x]);

Never use a name (in this case xi) and an indexed version (here xi[1] and xi[2] ) in the same worksheet.
Use other names for xi[1] and xi[2], e.g. xi1 and xi2. 
Try the following short session and you will see why:
 

restart;
xi := H[a]^2+(1+K)/K[p];
xi[1];
xi[1] := M*(1+K)*theta[1]^3/(K*xi)-(2+K)*theta[1]/xi;
xi;
eval(xi);

There is no similar problem for theta, since only indexed versions appear.
Next point: On my computer the integration of u doesn't finish before my patience runs out.
This integration is basically very simple. So start your worksheet like this (after having made the name changes mentioned).
 

restart;
h := x-> 1+x*tan(theta)/delta+phi*sin(2*Pi*x);
u := S*(xi1*cosh(theta[1]*y)*sinh(theta[2]*h(x))-xi2*cosh(theta[2]*y)*sinh(theta[1]*h(x))-L)/(L*xi)-1;
q := int(u, y = -h(x) .. h(x));

My last suggestion. Change the last two lines to

Px:=((q+2*h(x))*L*xi)/(F(x));
Deltap:=int(Px,x=0..1);

Copying your ode I find that you have a sequence of length 2 on the left (and a zero on the right).
Your ode is:

ode:=((0.9181604810e-6*cos(phi(t))^3+0.4213239281e-4*cos(phi(t)))*(diff(diff(phi(t), t), t))-0.9181604810e-6*(diff(phi(t), t))^2*cos(phi(t))^2*sin(phi(t))+0.2166831514e-7*(diff(phi(t), t))^2*sin(phi(t))^3*cos(phi(t))^2/(-0.2359970352e-1*cos(phi(t))^2+1)-0.5976753198e-5*(-.1536219500*(diff(diff(phi(t), t), t))*sin(phi(t))/(-0.2359970352e-1*cos(phi(t))^2+1)^(1/2)-.1536219500*(diff(phi(t), t))^2*cos(phi(t))/(-0.2359970352e-1*cos(phi(t))^2+1)^(1/2)+0.3625432474e-2*(diff(phi(t), t))^2*sin(phi(t))^2*cos(phi(t))/(-0.2359970352e-1*cos(phi(t))^2+1)^(3/2))*sin(phi(t))*(-0.2359970352e-1*cos(phi(t))^2+1)^(1/2)*cos(phi(t))+0.205427606e-5*(-.1536219500*(diff(diff(phi(t), t), t))*sin(phi(t))/(-0.2359970352e-1*cos(phi(t))^2+1)^(1/2)-.1536219500*(diff(phi(t), t))^2*cos(phi(t))/(-0.2359970352e-1*cos(phi(t))^2+1)^(1/2)+0.3625432474e-2*(diff(phi(t), t))^2*sin(phi(t))^2*cos(phi(t))/(-0.2359970352e-1*cos(phi(t))^2+1)^(3/2))*cos(phi(t))+0.5376768250e-3*(diff(phi(t), t))*cos(phi(t)), 

0.1836916543e-3*cos(phi(t))^2*(diff(diff(phi(t), t), t))-0.1836916543e-3*sin(phi(t))*cos(phi(t))*(diff(phi(t), t))^2-0.2821907012e-4*(diff(phi(t), t))^2*sin(phi(t))^2*cos(phi(t))/(-0.2359970352e-1*cos(phi(t))^2+1)^(1/2)+0.7783642454e-2*(-.1536219500*(diff(diff(phi(t), t), t))*sin(phi(t))/(-0.2359970352e-1*cos(phi(t))^2+1)^(1/2)-.1536219500*(diff(phi(t), t))^2*cos(phi(t))/(-0.2359970352e-1*cos(phi(t))^2+1)^(1/2)+0.3625432474e-2*(diff(phi(t), t))^2*sin(phi(t))^2*cos(phi(t))/(-0.2359970352e-1*cos(phi(t))^2+1)^(3/2))*(-0.2359970352e-1*cos(phi(t))^2+1)+0.1210411733e-2*(diff(diff(phi(t), t), t))*sin(phi(t))/(-0.2359970352e-1*cos(phi(t))^2+1)^(1/2)+0.1210411733e-2*(diff(phi(t), t))^2*cos(phi(t))/(-0.2359970352e-1*cos(phi(t))^2+1)^(1/2)-0.2856535802e-4*(diff(phi(t), t))^2*sin(phi(t))^2*cos(phi(t))/(-0.2359970352e-1*cos(phi(t))^2+1)^(3/2)+1.553675255*(50*exp(-2.601456816*(phi(t)-4*Pi*trunc((1/4)*(phi(t)+Pi)/Pi)-.36)^2)+2)*Pi*cos(phi(t))) = 0;

Surely that comma ougth to be something else or you only want one of  the two?
The error comes in at the definition of Bewegungsgl:
nops([lhs(Bewegungsgl)]); # answer 2
In fact already here:
nops(GleichungKräfte[2]);

You can do as follows:

restart;
eq:=diff(x(t), t$2) = .8/(x(t)^3*exp(0.02*int(1/x(s), s = 0 .. t-tau)))-1/x(t)^2;
ics := x(0) = 1, D(x)(0) = 0;
eq_y:=y(t)=int(1/x(s), s = 0 .. t-tau);
ic_y:=eval(eq_y,{t=0,x=1});
sys:={diff(eq_y,t),subs(rhs(eq_y)=lhs(eq_y),eq)};
res:=dsolve(eval(sys union {ics,ic_y},tau=0.1),numeric);
plots:-odeplot(res,[t,x(t)],0..30);

After giving your system the name sys I did:
 

g:=indets(sys)=~1;
fsolve(sys,g);

and quickly got
{c[0][0] = 12.89605073, c[0][1] = -6.470667865, c[0][2] = 1.752111454, c[0][3] = -.1774943190, c[1][0] = 9.737598780, c[1][1] = -2.590154874, c[1][2] = .3402465368, c[1][3] = -0.1700267263e-1, c[2][0] = 7.032179265, c[2][1] = -1.179016309, c[2][2] = 0.9696822774e-1, c[2][3] = -0.3121342845e-2}
which agrees well with your second real solution.
The title of your question says that fsolve gives an error with a range. It didn't in your worksheet; it just returned unevaluated, i.e. it didn't succeed in finding a solution.
Often I find that using a guess instead of a range (or ranges) is surer of success.

Even the version in the help page for semitorus in Maple 2016.2 looks weird.
Using capped=false gives you a nice plot, but with no cap.
Here is what I got in Maple 2015.2 using your command as given, i.e.

plots:-display(plottools:-semitorus([0, 0, 0], 0 .. Pi, 1, 2), lightmodel = light4, orientation = [-140, 60], scaling = constrained, style = patchnogrid);

If you turn the plot (in Maple, not here) you will see that it looks like a donut cut in half horizontally and with a piece of paper at the bottom. The "paper" is due to the default capped=true.

I get the same thing in Maple 18 and 15, while in Maple 12 I get something quite different:

Obviously here the angle means something different than in the other versions.

To get something like this Maple 12 version in Maple 2016.2 you could try tubeplot (in the plots package. Here I take 3/4 of the donut cut vertically, but not capped:

tubeplot([cos(t),sin(t),0],t=0..3*Pi/2,radius=1/2,numpoints=60,tubepoints=20,scaling=constrained);

I just checked in Maple 8. It works as Maple 12.
I wonder if there was a documented change in Maple 13 or 14 about the meaning of the angular range (the second argument) in semitorus.
My suspicion is that this is a bug introduced in Maple 13, 14, or 15.
I notice that the help pages are the same.
In Maple 2016 you find
"The semitorus command creates a three-dimensional plot data object, which when displayed is a section of the torus centered at [x,y,z], whose meridian has radius r and distance from the center of a meridian to the center of the semitorus is R.  The semitorus starts at radian angle a and ends at radian angle b."
In Maple 12 you find:
"The semitorus command creates a three-dimensional plot data object, which when displayed is a section of the torus centered at [x,y,z], whose meridian has radius r and distance from the center of a meridian to the center of the semitorus is R. The semitorus starts at radian angle a and ends at radian angle b."

That looks the same to me.
I will submit an SCR.