I tried the model suggested by Carl Love and modified by Kitonum in NonlinearFit in the Statistics package.
 

restart;
M:=ImportMatrix("C:/Users/Bruger/Downloads/Book1.xlsx",datatype=float[8]);
M[1..10,..];#Having a look
mx:=max(M[..,1]); #maximal x-value
F:=Statistics:-NonlinearFit(A*x^n*sin(B*x),M,x,initialvalues=[A=0.015,B=280,n=1.5],output=solutionmodule);
F:-Results();
res:=F:-Results("leastsquaresfunction");
sq1:=F:-Results("residualsumofsquares"); #For comparison with Kitonum's function
f:=x->0.015*x^1.5*sin(280*x); #Kitonum
Y:=f~(M[..,1]);
sq2:=add((Y[i]-M[i,2])^2,i=1..numelems(Y)); #"residualsumofsquares" for f.
p1:=plot(M,style=point):
p2:=plot(res,x=0..mx,color=blue,numpoints=30000):
plots:-display(p1,p2,size=[1800,default]);
p3:=plot(f(x), x=0..mx,color=blue,numpoints=30000):
plots:-display(p1,p3,size=[1800,default]);
plots:-display(p1,p3,size=[1800,default],view=[mx-.1..mx,default]);
plots:-display(p1,p2,size=[1800,default],view=[mx-.1..mx,default]);
################# Introducing 2 extra parameters:
F2:=Statistics:-NonlinearFit(A*x^n*sin(B*x^m+phi),M,x,initialvalues=[A=0.015,B=280,n=1.5,m=1,phi=0],output=solutionmodule);
res2:=F2:-Results("leastsquaresfunction");
F2:-Results("residualsumofsquares"); #0.0008
p22:=plot(res,x=0..mx,color=blue,numpoints=30000):
plots:-display(p1,p22,size=[1800,default]);

The result 'res' was

Kitonum's function plot in some way looks better by hitting the tops, but the residual sum of squares is 10 times larger than that from NonlinearFit. Also the closeups at the end:

plots:-display(p1,p2,size=[800,default],view=[mx-.1..mx,default],caption="NonlinearFit");
plots:-display(p1,p3,size=[800,default],view=[mx-.1..mx,default],caption="Kitonum's function");

are revealing:

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.