1. In order to compute such integral, Maple would need to know the Frullani integral, and it does not. Not even the classical (usual) one found in any calculus textbook., e.g.

int( (h(a*x) - h(b*x))/x, x=0..infinity )
where
h := x -> sin(sin(x))/x;
This should not be difficult to implement.

2. In your example one needs a generalized Frullani theorem (less known), because the limit of h at infinity does not exist. Namely, instead of limit(h(x), h=infinity), it is enough to exist
limit( int(h(x),x=x0..X)/X, X=infinity) for some x0>0. Actually in your example h is periodic, and this limit is int(h, 0..P)/P, P being the period (P = Pi in this case).

3. To compute numerically such integrals Maple would need good numerical methods for oscillatory integrals.

4. Please note that you used wrong formulae for computing A and B. The fact that the results are correct is just by chance: a-b=1. 
The correct A, B are:
J := int(abs(sin(x)), x = 1 .. T) assuming T>4:
B := limit(J/T, T = infinity);
K := x*int(abs(sin(t))/t^2, t = x .. 1) assuming x>0,x<1;
A := limit(K, x = 0, right); # or simply  A = abs(sin(0)) = 0 in this case.

You could use Statistics:-Fit to find the function. If you know the theoretical model, use it; otherwise you may try a polynomial.
This should be enough.
If you really want to remove some points, the best candidates are those with large residuals.

mseries:=proc(f,X,n)
local t;
eval( f, `=`~(lhs~(X), rhs~(X)+t*~(lhs~(X)-rhs~(X))) );
eval(convert(series(%, t, n), polynom),t=1);
end;

# test
mtaylor(sin(x^2+y^2), [x=1,y=0], 3) = mseries(sin(x^2+y^2), [x=1,y=0], 3):
simplify(%);

Note that y=0 cannot be shortened to y  in this version.

Is you use numerical integration, i.e.

u := proc (x, t) options operator, arrow;
   add(int(sin(k*lambda*x)*g, x = 0 .. 1,numeric))*exp(-k*lambda*t)*sin(k*lambda*x), k = 1 .. n)  end proc;

the code (without threads) is 40 times faster.

The first bug is that Mode returns the set of critical points (actually float approximations) instead of the maximum point.
Trying to isolate the error, a bug in solve related to LambertW function appeared. Here is the commented Maple code.

with(Statistics):
X := RandomVariable(Normal(-2,1)):
Y := RandomVariable(Normal(2,1)):
r    := 4/10:      # r=1/2 is similar
f:= (1-r)*PDF(X,t)+r*PDF(Y,t);

f__Z := unapply(f, t);
Z    := Distribution(PDF=f__Z):
Mode(Z);
   Warning, solutions may have been lost 

maximize(f,location);  #does not work, returns unevaluated!
maximize(f ,t=-10..10,location);  # wrong answer,  the max is not at t= -2

f1:=diff(f,t):
solve(f1,t);   # WRONG!!!  returns the result of Mode
s:=allvalues(%);
eval(f1, t=s[1]);  evalf(%);  #  so, not a root!
     

        -0.2141283611e-3

ans:=fsolve(f1,t=-3..-1);
      ans := -1.999102417
evalf(eval(f1, t=%));
     

## So, the true Mode is ans  [as float approx].

In Maple only the first question has to be answered, the rest are automatic (almost).
See:
?sum   ?exp   ?infinity

Hint: start with
assume(delta <= 0);
(for convergence).

Q:=Quantile(X, z);
simplifies to
RootOf(_Z^5-30*_Z^3-875*_Z+10000*z-5000)

But this is not quite correct; the polynomial in _Z may have more than one real root and the correct one is in -5..5). The plot in the first case is for another branch. In the second case, Q is not evaluated and the plot is correct but very slow (because probably fsolve is called)!

The fast approach is to choose the correct branch in RootOf.
RootOf(op(Q), -5..5);
plot(%, z=0..1);
 

mtaylor(sin(x+y)+y, [x=Pi, y=0], 6);

sinc:= x -> sin(x)/x:
sinc(0):=1: sinc(0.):=1.:
n:=10;  #1000 if you insist
a:=5.0; h:=a/n;
Z:=Matrix(2*n+1, (i,j) -> sinc((i-1-n)*h)^2*sinc((j-1-n)*h)^2, datatype=float[8]);

# note that sinc is even, so Z is symmetric.

You can convert using IrfanView (free)
http://www.irfanview.com/main_formats.htm

Explore(plot([(-y0-x0)*exp(-2*t)+(y0+2*x0)*exp(-t),
               -2*(-y0-x0)*exp(-2*t)-(y0+2*x0)*exp(-t), t=0..5]),
        x0=1..5., y0=1..5.);

evalb(irem(390,7)=0);
   false

Q:=[X,Y]:
Ps:=proc(f,g)
local k;
global P,Q; # not really necessary
add(diff(f,P[k])*diff(g,Q[k]) - diff(f,Q[k])*diff(g,P[k]),k=1..nops(Q))
end: 

Ps(P[1],X) ;

with(plots):with(plottools):
X:=cos(t):Y:=sin(t):
vol:=Pi*X^2*(2*Y) + 2*Pi*X^2*(1-Y)/3:
sol:=maximize(vol,t=0..Pi/2,location):
t:=eval(t, op([1,1],sol[2]));#evalf(%)*180/Pi;
p1 := display(polygon([[[-X, Y], [0, 1], [X, Y], [X, -Y], [0, -1], [-X, -Y], [-X, Y]]]), color = red, thickness = 2):
p2 := implicitplot(x^2+y^2 = 1, x = -1 .. 1, y = -1 .. 1):
display(p1, p2);
simplify([X,Y]);

#[X,Y]

It is not fsolve, it is your system which contains

[c0 = 0, ... , (...)/c0 + ... = 0, ...]

Of course it cannot be solved. Replace ":" by ";" to see what happens.

PS. Why do you use such strange constructs?
F(eta) := ...
...
F(eta) := unapply( F(eta), eta);
Sintactically it is correct but now F is an empty procedure and in its remember table F(eta) is also a procedure. It works, but why so?