Any 3D vector graphics is terrible or simply not working. The old Classic Interface is a little better but this interface is now absent for a 64 bit OS.
I have posted in the past some examples here but I cannot find them.
It seems that Maple considers this (very old!) problem as irrelevant. I use Maxima (free) for vector 3d graphics.
Here is how a simple plot
plot3d(x^2+y^2, x=-1..1,y=-1..1);
is exported as encapsulated postscript (in pdf it is exported as bitmap!):

P:=(u,v)->[u,v,-(u^2+v^2)]:
r:=(u,v,t)-> `if`(-2*u*cos(t)-sin(t)<0,0.1,1):
g:=(u,v,t)-> `if`(-2*u*cos(t)-sin(t)<0,0.1,0):
b:=g:
plots:-display(
seq( plot3d( P(u,v), u=-1..1, v=-1..1,  color=[r(u,v,t),g(u,v,t),b(u,v,t)]), t=0..2.*Pi,Pi/20.),
insequence=true);

Or, you can use Explore:

Explore( plot3d( P(u,v), u=-1..1, v=-1..1,  color=[r(u,v,t),g(u,v,t),b(u,v,t)]),  t=0..2*Pi, animate, loop,autorun) ;

If d(x,y,z) is the divergence, the integral is:

int( d(x,y,z), [z=x^2+y^2..-2*x, y=-sqrt(-2*x-x^2)..sqrt(-2*x-x^2), x=-2..0]);

Plot:
plot3d(  [x^2+y^2,  -2*x],  x=-2..0, y=-sqrt(-2*x-x^2)..sqrt(-2*x-x^2));

plot3d of
s:=[sin(phi)*cos(theta), sin(phi)*sin(theta), cos(phi)];
for phi, theta satisfying

{op(s <= ~ a), op(s >= ~ -a),  phi<=Pi, Phi>=0, theta<=2*Pi, theta>=0}

Unfortunately solve is not able to solve this system of inequalities, so you will need to do it by hand (which is not difficult; you may also use the symmetry).

The fact that
asympt(Zeta(1/2+I*y),y,6): simplify(%);
     1+2^(-1-(2*I)*y)+3^(-1/2-I*y)+2^(-1/2-I*y)

"works"  is just a mistake.  Please note the absence of the O(...) term.
Note also that
MultiSeries:-asympt(Zeta(1/2+I*y),y,6);
==> Error, (in MultiSeries:-multiseries) unable to compute multiseries

N:=1000:
f:=n ->evalf(frac(n*sqrt(2))):
g:=n ->evalf(frac(ln(n))):
X:=Vector(N,i->i):
Yf:=Vector(N,f):  Yg:=Vector(N,g):
plots:-display( Array([
plot(X,Yf, style=point), plot(X,Yg, style=point)]) );

Here is a crude version of a Frullani proc.

Fru:=proc(f,x,a,b,x0:=1)
local A,B,X,JA,JB;
B:=limit(f, x=infinity);
if type(B,`..`) or has(B,{limit,undefined}) then
  JB:=int(f,x=x0..X) assuming X>1000*(x0+1);
  B:=limit( JB/X, X=infinity)  fi;
A:=limit(f,x=0,right);
if type(A,`..`) or has(A,{limit,undefined}) then
  JA:=int(f/x^2,x=X..x0) assuming X>0,X<x0;
  A:= limit(X*JA, X=0, right) fi;
Int( (eval(f,x=a*x) - eval(f,x=b*x))/x, x=0..infinity) = (B-A)*ln(a/b)
end;

#tests

Fru( abs(sin(x)),x,2,1);
Fru( sin(sin(x))/x, x, a,b);

Edit: forgot an "x" in Int

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