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

It will be impossible to represent such a huge object.

n:=30:
A:=simplify(LinearAlgebra:-RandomMatrix(n)+x):
B:=A^(-1):
length(B);

    2334559

If you have enough time, try then  n:=100.

restart;
f := z -> MeijerG([[], [1]], [[0, 2], []], z):
e :=      MeijerG([[], [1]], [[0, 2], []], z):
plots:-display(plot(f, -2.5..2,color=red), plot(e, z=-2.5..2,color=blue));

There are infinitely many solutions.

solve(z^(1+I) = 1, z, allsolutions):
evalc(%);

So, you don't need evalf.

BTW, seq(print(...),...)  is strange (and the output is wrong; what version are you using?). Why not:

interface(rtablesize=infinity):
Vector([seq(evalf(evalf[d](u)), d = 10 .. 3010, 300)]);

I think that readstat being  an I/O command, it may need time to complete before calling it again.
Inserting some Sleeps makes it work.
dt depends on the speed of the computer.

dt:=0.5:
fun:=proc(n)  local a, b, i, x;
a:=readstat("insert a");
Threads:-Sleep(dt);
b:=readstat("insert b");
for i to n do
Threads:-Sleep(dt);
x[i]:=readstat("insert x");
od;
end proc:
fun(3);

sqrt (the principal branch) is af course well defined, the branch cut is arbitrary (-oo,0); it could be any "arc" from 0 to oo.

But it is possible to define a custom sqrt. Here is one with a branch cut along negative imaginary axis.

restart;
mysqrt:=z -> piecewise(argument(z)>-Pi/2,sqrt(z),-sqrt(z)):
 

limit(mysqrt(-1-t*I), t=0, left), limit(mysqrt(-1-t*I), t=0, right); #continuous
                              I, I
limit(sqrt(-1-t*I), t=0, left), limit(sqrt(-1-t*I), t=0, right);     #discontinuous
                             I, -I
limit(mysqrt(-I-t), t=0, left), limit(mysqrt(-I-t), t=0, right);     #discontinuous

   
limit(sqrt(-I-t), t=0, left), limit(sqrt(-I-t), t=0, right);        #continuous        

(1/2)*sqrt(2)-I*sqrt(2)*(1/2), (1/2)*sqrt(2)-I*sqrt(2)*(1/2)

simplify(convert(z, arcsin)) assuming x>-a,x<a;

works in Maple 2017.

Download loss.mw