With a and b as in your worksheet the simplest way seems to be:

Equate(a.b,<(0$4)>);
res:=solve(%);
eval(a.b,res); #Just a check

But you probably meant to require that the unknowns are real. Correct?

If so then you will have 8 linear equations with 7 unknowns. That system is likely not to have any solution. Indeed it appears so:

eqs:=Equate(a.b,<(0$4)>);
eqs1:=evalc(Re~(eqs));
eqs2:=evalc(Im~(eqs));
solve({op(eqs1),op(eqs2)});

Since your numbers are decimal numbers (floats) the existence of a solution to your overdetermined system will be sensitive to roundoff in input etc.

Thus you may want to settle for the best substitute for a solution in the least squares sense:

nms:=convert(indets(b),list);
A,B:=LinearAlgebra:-GenerateMatrix([op(eqs1),op(eqs2)],nms);
X:=LinearAlgebra:-LeastSquares(A,B);
LinearAlgebra:-Norm( A . X - B ,2); # good if small
eval(a.b,nms=~convert(X,list)); #Another check

You can solve numerically.
The heat equation immediately smoothes the initially given function and your boundary conditions bc:=f(t,0)=0,f(t,1)=1 together with the initial condition ic:=f(0,x)=piecewise(x=0,1,0) disagree at (t,x) = (0,0) and (t,x) = (0,1). Thus you can expect to see problems initially in the numerical solution.
I wonder, however, if you didn't intend ic to be ic:=f(0,x)=piecewise(x=1,1,0) instead? That wouldn't help for an exact solution, but you won't see any numerical problems:
 

sol:=pdsolve(pde,{bc,ic},numeric,spacestep=0.05,timestep=0.001);
sol:-animate(t=0.4,frames=100);

Use numerical solution:
 

restart;
ode:=diff(y(x),x,x)=5*exp(-10/diff(y(x),x)); 
bcs:=y(0)=0,y(15)=2;
res:=dsolve({ode,bcs},numeric);
plots:-odeplot(res,[x,y(x)]); # Looks like a straight line y=2*x/15
plots:-odeplot(res,[x,y(x)-2*x/15]); # basically 0
odetest(y(x)=2*x/15,ode); # Almost zero, but not quite

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