After studying Maple such things will be possible. E.g.

funny := module() option package;  export `<`;
   `<`:= proc(a::list(realcons), b::realcons) option overload;
              select(t -> t<b, a);
         end proc;
end module:

with(funny):
a:=[1,2,3,4,5]:
b:=a<4;
                         b := [1, 2, 3]

Suppose you want F[i, j,...](x, y,...)  to return  x^i + y^j + ...

F:=proc()
  local ind:=op(procname); 
  add( [args] ^~ [ind] )
end:

F[i,j,99](x,y,z);
    x^i+y^j+z^99

Of course some parameter checks in the procedure are recommended.

b:=A[..,6]:
C:=A[..,[1..5,7..11]]:
(C^(-1).b)[6];

(simplify@expand@convert)( abs(u+v)^2+abs(u-v)^2, conjugate );

1.   As documented.

2. Because is tries to be smart.
Also:
is(0.5 = 1/2);
     true

3. Because is is not perfect.
Also:
is(exp(1)+Pi, rational);
    false
Should be FAIL, the result is not known (yet)!

f:=exp(t^(1/3)-2*t^(7/4)):
s:=series(f,t,6);

sort(s,t,ascending);

The timing seems to be similar.

t:=time[real]():
L:=combinat:-permute([1,1,1,1,0,0]):
C:=Iterator:-CartesianProduct([seq(1..15)]$6):
n:=0:
T:=table():
for c in C do
if add(L[c[i]][1],i=1..6)<>4 or
   add(L[c[i]][2],i=1..6)<>4 or
   add(L[c[i]][3],i=1..6)<>4 or
   add(L[c[i]][4],i=1..6)<>4 or
   add(L[c[i]][5],i=1..6)<>4 or
   add(L[c[i]][6],i=1..6)<>4 then next fi;
n:=n+1;
T[n]:=Matrix([seq(L[c[i]],i=1..6)]);
od:
'num'=n, 'time'=time[real]()-t;

        num = 67950, time = 65.329

Edit. We gain about 2 seconds if one of the 6 tests inside if is removed.

@shimaa sadk 

The inner sum can be computed symbolically. Do it.
Fot the outer one use evalf(Sum(...,T=0..infinity);

Unfortunately evalf for GAMMA and also Zeta is VERY slow for large arguments. E.g.

evalf( GAMMA(1+10^20*I));

The evalf operations are known to be hard/impossible to interrupt.

It seems that for infinite complex limits Maple applies some nonvalid transforms (fractional powers).

K:=Int(1/(s^2+1), s = 1-I*infinity .. 1+I*infinity):
eval(K, infinity=1000):
evalf(%);
    -9.*10^(-13)+0.1999998661e-2*I

Transforming by hand the complex limits also works (IntegrationTools:-Change  fails):

J:=Int(  1/((1+I*t)^2+1)*I, t=-infinity..infinity):
evalf(%);
    1.110223025*10^(-16)*I

[1.] p1 is not function (in mathematical sense, i.e. a procedure in Maple). It is an expression.

[2.] diff works, but not as you expect. You should be aware that p1 could be nowhere differentiable; this depends of a(t).
For example if a(t) = 1 for a rational t and 0 otherwise then you cannot differentiate p1.

[3.] p1(1) is a nonsense, see [1.].  Think e.g. at the meaning of  sin(t)(1).

[4.] The simplification is indeed strange, the first vector has disappeared.

[5.] See again [1.]

int~(int~(a, t=0..t), t=0..t) + v0*t + r0;

or even shorter:

int~(a, t=0..t, t=0..t) + v0*t + r0;

Actually the system is linear, so it can be solved exactly even without dsolve.
The solutions are huge, so must be ploted only near 0.

Eigenvalues uses special algorithms to obtain best accuracy for the result (float[8] entries).
Computing the inverse of RS introduces large roundoff errors.
Note also that RS is almost singular; its determinant is ~ 10^(-1208);