Not answering your actual question just giving the value assuming that the parameters are real and different from zero in which case we may safely assume that they are positive.
A:=Int(exp(-xi^2*b^2)*cos(xi^2*a*t)*cos(xi*x), xi = 0 .. infinity);
value(A) assuming positive;

@ESchiltkamp I'm now very certain that the exact result returned by dsolve for your ode named 'somFy' is wrong. I suspected that, but am now sure. The difference in your case between the correct solution and the dsolve-solution is small enough to make the result look reasonable. However, it is wrong.

It helps to look at a version with simpler constants:
restart;
ode:=-1-diff(y(t), t)^2*signum(diff(y(t), t)) =diff(y(t), t, t);
ics:=y(0) = 0, D(y)(0) = 1;
res:=dsolve({ode,ics});
rr:=odetest(res,ode); #Ought to be zero
plot(rr,t=-Pi/4..3*Pi/4,-.1..1,thickness=3);
#At t = Pi/4 the derivative of res changes sign since cot changes sign at Pi/2. After that rr is nonzero.
diff(res,t);
#The result res is the same as
u:=ln((sin(t)+cos(t))^2)/2;
#thus equal to ln(abs(sin(t)+cos(t))) for t real.
#Now ignoring signum:
ode2:=eval(ode,signum=1);
res2:=dsolve({ode2,ics});
odetest(res2,ode2);
##########
So res and res2 agree when sin(t)+cos(t) > 0. However, that is the case for t=0. Since sin(t)+cos(t) > 0 on -Pi/4 .. 3*Pi/4 and  tends to 0 as t tends to 3*Pi/4 from the left and also as t tends to -Pi/4 from the right the solutions both have maximal domains equal to the open interval -Pi/4 .. 3*Pi/4. That res makes sense for t > 3*Pi/4 is irrelevant because of the singularity at 3*Pi/4.
##########
#Solving ode numerically:
resN:=dsolve({ode,ics},numeric);
evalf(Pi/4);
plots:-odeplot(resN,[t,y(t)],-Pi/4..5); pnum:=%:
##########
#By solving ode2 up to t=Pi/4 (where y ' (t) = 0) and continuing with
ode3:=eval(ode,signum=-1);
#with the approbriate ics at t=Pi/4 it is starightforward to finf the actual exact solution:
RES := y(t) = piecewise(t < Pi/4, ln(sin(t)+cos(t)), t+Pi/4+(3/2)*ln(2)-ln(exp(2*t)+exp((1/2)*Pi)));
##Plotting
plot(rhs(RES),t=-Pi/4..5); pex:=%:
##Displaying RES with the numerical solution
plots:-display(pnum,pex); #They fall on top of each other.

I shall file a bug report (SCR).

@ESchiltkamp Changing alpha to 20 degrees does give us an answer to displacementy which doesn't contain a RootOf. However, I believe it is wrong. First of all it is suspicious that there is no cutoff like what happens in a piecewise. Secondly, odetest gives an identically zero result from t=0 to some t0, after that it continues continuously to nonzero values (untill it blows up, which is not the problem).
displacementy := dsolve({ics2, somFy});
rrr:=odetest(displacementy,somFy): #Testing, rrr ought to be zero
plot(eval(rrr,v0=100),t=0..0.3,0..1);
###############
I tried a different approach. Purely numerical.
Digits:=15:
res:=dsolve(sys,numeric,parameters=[v0],output=listprocedure);
X,Y:=op(subs(res,[sx(t),sy(t)]));
res(parameters=[5]); #Setting the v0 parameter for a test
odeplot(res,[sx(t),sy(t)],0..10,numpoints=100);
odeplot(res,[[t,sx(t)],[t,sy(t)]],0..3,numpoints=100);
odeplot(res,[[t,diff(sx(t),t)],[t,diff(sy(t),t)]],0..3,numpoints=100);
#Now setting up a way to find the speed v0 necessary.
p:=proc(T,v0) option remember; res(parameters=[v0]); [X(T)-distance,Y(T)] end proc;
px:=proc(T,v0) p(_passed)[1] end proc;
py:=proc(T,v0) p(_passed)[2] end proc;
p(1,5); #Testing
sol:=fsolve([px,py],[3,3e12]); #A huge speed guess
res(parameters=sol[2..2]); #For safety I set the parameter v0
odeplot(res,[sx(t),sy(t)],0..sol[1]);
X(sol[1]); #This value ought to be 2.5, but it is somewhat close
#I suspect that there are some numerical problems since the input of [3,3e12] into fsolve returns a listwith the v0-guess unchanged. It is probably due to the huge value of v0.
Trying with
sol:=fsolve([px,py],[3,3.18e12]);
results in a value closer to distance=2.5.
Question: Are the parameters rho and Cd in your problem correct?
You have an enormous resistance requiring an enormous starting speed for this very light ball to hit the ground at 2.5 m.
#################
After seeing your corrected version I tried changing the procedure p above to
p:=proc(T,v0) option remember;
  local L,eps;
  eps:=10^(-Digits);
  res(parameters=[v0]);
  L:=[X(T)-distance,Y(T)];
  if L[1]^2+L[2]^2<eps then print([T,v0],L) end if;
  L
end proc;
px:=proc(T,v0) p(_passed)[1] end proc;
py:=proc(T,v0) p(_passed)[2] end proc;
##That way I get away from the somewhat unwieldy fsolve. I force it to print out results that are good enough for me, but not not good enough for fsolve.
sol:=fsolve([px,py],[0..1,0..7]);

@Carl Love The statement "Anything inside the quotes is treated as a variable, regardless of any other meaning(s) that the characters have." should be modified in view of the results of
`sin`(Pi);
sin(`Pi`);
`sin`:=5; #You won't be allowed to do this
and similar.
The help page ?names says it like this:
"Any valid Maple name formed without using left single quotes is precisely the same as the name formed by surrounding the name with left single quotes. "
sin and Pi are (also) names.

@shiMusa
Yes, the point is not that the numbers should be rationals. The point is that they should not be floats.
The following is straightforward:
restart;
M := Matrix(2, 2, {(1, 1) = a+2.5*I, (1, 2) = 1-I*a, (2, 1) = 4, (2, 2) = a});
M3:=eval(M,a=0.5+0.5*I);
convert(M3,rational);
##
This version is more subtle:
restart;
M := Matrix(2, 2, {(1, 1) = a+2.5*I, (1, 2) = 1-I*a, (2, 1) = 4, (2, 2) = a});
a:=0.5+0.5*I;
M;
convert(M,rational); #No effect
#However, elementwise application of convert works fine:
convert~(M,rational); #Notice the tilde, which implies elementwise application of convert
#To get an idea of what is going on, try debugging:
debug(`convert/rational`);
convert(M,rational);
#Notice the input and the contents of 'flts'.
# Continue with
convert~(M,rational);
#Notice the difference!
#I guess the conclusion is that you should use elementwise conversion.
# To get back to the normal state use undebug:
undebug(`convert/rational`);
##Now try at the element level: first row, first column:
convert(M[1,1],rational);
# We can continue with
indets(M,name);
%;
indets(M,float);

@oldstudent Actually I edited your expression eval( x^3-3*x^2-24*x+8, [x = -2])  to what I wrote in my answer above. Notice in particular the tilde (~) which is used twice in this example. 
The first that happens in
eval~( x^3-3*x^2-24*x+8, x=~[-2,4,6,9]);
is that x=~[-2,4,6,9] becomes [x = -2, x = 4, x = 6, x = 9], which you can try like this:
x=~[-2,4,6,9];
This is use of elementwise `=`.
After that eval is applied elementwise to the list of x-equations as in
eval~(x^3-3*x^2-24*x+8,[x = -2, x = 4, x = 6, x = 9]);

I don't understand what you want to do.
Your worksheet contains the assignment
E := sum((V[i]-x[i]*((1+b)/(b*x[i]^4+1))^((1+b)/(4*b)))^2, i = 1 .. 78);
and nothing else.
What do you want to do with it?

Apparently the singularity is a problem.
Try
plots[contourplot](1/(x^2+y^2), x=-1..1, y=-1..1,contours=[seq(1..50,5)],grid=[100,100]);
plots[contourplot3d](1/(x^2+y^2), x=-1..1, y=-1..1,view=0..100);

You may try
?2DMathshortcutkeys

@Carl Love It appears that emmantop wants to (has to) write his own algorithm:

http://www.mapleprimes.com/questions/202416-How-Do-I-Write-Maple-Programme-For-Nonlinear-Ode

In that link the tag 'homework' appears. So I guess it is 'has to'.

Have you already got a program for the linear case? There the set of equations for the discrete version can be solved by LinearSolve immediately.
That is not so in the nonlinear case, where you will need to use Newton's method, but LinearSolve can still be used in each iteration.
The initial work in making a procedure handling the nonlinear case is basically the same as for the linear case. So making a procedure handling the linear case first makes sense to me.
I would make the procedure so it takes a system of first order as input.
 

@J4James In your test case this could be done like this:
restart;
eq1:=diff(f(y), y$4)-(diff(f(y), y$2));
bcs:=f(h1) = (1/2), f(h2) = -1/2, D(f)(h1) = -1, D(f)(h2) = -1:
h1:= 1+cos(x):
h2:=-1-cos(x+g):
db:=eq1,bcs:
d1 := subs(g=1,[db]);
F:=proc(xx,yy) local fp; global x,y;
   fp:=subs(dsolve(eval(d1,x=xx), numeric,output=listprocedure),f(y));
   fp(yy)
end proc;
   
F(0,.3); #Test of f(x,y)
plot3d(F,0..1,-2..2,grid=[20,20]);

##################
While the method above works it is unnecessarily slow: In doing the plot there will be n^2 calls to dsolve if the grid is [n,n]. The number of calls can be reduced to n if we change to:
Fx:=proc(xx) option remember; global y,x;
   if not type(xx,numeric) then return 'procname(_passed)' end if;
   subs(dsolve(eval(d1,x=xx), numeric,output=listprocedure),f(y))
end proc;
Fx(0)(.3); #Test
plot3d(Fx(x)(y),y=eval(h2..h1,g=1),x=-2..2);

 

@k20057 5
1. The Maple package Statistics has a procedure called NonlinearFit. I'm using the long name Statistics:-NonlinearFit for NonlinearFit. An alternative is to issue the command with(Statistics); after which you can use the name NonlinearFit instead.
2. NonlinearFit handles models like f(E)=B/(A*exp(E/kT)+1) containing unknown parameters, here B,A,kT. E is the variable. M is a matrix containing the experimental data with the E's in the first column and the values of the corresponding experimental f's in the second, here I shall refer to those as fexp(E).
3. NonlinearFit tries to find values for A,B,kT which minimizes the sum of squares of f(E)-fexp(E) over the given E's, i.e. add( (f(E[i])-fexp(E[i]) )^2, i=1..N) is attempted minimized.
4. The input to NonlinearFit is the sequence: Model (f(E) ), data (M), variable (E).

Edited: My remark about the value of A was due to a stupid typo.

Looking in a book (Kittel) on Solid State Physics I found a version of the distribution written like this:
f(E)=1/(exp(E-E_F)/(k*T)+1).
Written like that A would be given by A=exp(-E_F/(k*T)), thus positive.
So I wonder about the images you posted.

Maybe a normalization is needed (reflected in the parameter B):
f:=E->B/(A*exp(E/kT)+1);
Now you could either write down 3 equations using 3 points and solve for A,B,kT or (better) use all the points thereby getting an f that is best possible in some sense.
f:=E->B/(A*exp(E/kT)+1);
eq1:=f(0)=1.8;
eq2:=f(1)=1.6;
eq3:=f(2)=1.2;
solve({eq1,eq2,eq3},{A,kT,B});
plot(eval(f(E),%),E=0..9);
#The latter approach:
M:=Matrix([[0,1.8],[1,1.6],[2,1.2],[3,0.8],[4,0.4],[5,0.2],[6,0],[7,0],[8,0],[9,0]]);
res:=Statistics:-NonlinearFit(f(E),M,E);
plot([res,M],E=0..9,style=[line,point],symbolsize=25);

@blink23 Maybe ModuleLoad is what you need.
?ModuleLoad
So remove the protect command, make ModuleLoad a local and add these lines:
ModuleLoad:=proc() protect(i,j,k,epsilon); printf("Protecting i,j,k,epsilon\n") end proc;
ModuleLoad();

Will that be OK?