You must use parallel substitutions in ex:
ex := simplify(subs([ X = cos(theta)*X-sin(theta)*Y, Y = sin(theta)*X+cos(theta)*Y ], eq)); 

(You have defined the procedure f. It is then easier to use f(X-9/7, Y+8/7)  and similar for ex, instead of subs) .

It is interesting that for diff, the first argument can be almost anything, including a mathematical nonsense. For sets, lists, rtables it acts as expected (elementwise), but even strings are accepted. 
Note that series may produce errors.

f:=sin(x) + x^3*"a";
g:=sin(x) + [cos(x), x^3];
h:=sin(x) + [cos(x)*"a"] + {x^2*"b"^3 + "c"};

You can reduce the time to 0 because if a,b,... are >0 (as you have tried) then g must be a square.

The problem is generated by simplify.

(1 + (203808*exp(-(342569*t)/506))/131537)^(131537/203808);

simplify(%, sqrt);  # or simplify(%)

For example, with a smaller exponent:

simplify((1 + (203808*exp(-(342569*t)/506))/131537)^(3/7), sqrt); 

You should try to understand and explain what happens in the next lines:

restart;
P:=proc(a,b,c)
local S:=LinearAlgebra:-RandomMatrix(3, generator=rand(-10..10),shape=skewsymmetric),
      V0:=LinearAlgebra:-RandomVector(3, generator=rand(-10..10))^+,
      A:=(S-1)^(-1) . (S+1) . <a,0,0;0,b,0;0,0,c>;
A:=ilcm(entries(denom~(A),nolist))*A^+;
<V0, V0+A[1],V0+A[2],V0+A[3]>
end:
seq(P(3,5,14), 1..10);

Only a few miliseconds are needed for the 128 solutions.

restart;
ti:=time[real]():
XY:=seq(seq(`if`(irem(x*y,x+y)=0,[x,y],NULL), y=2..x-1),x=2..30):
N:=0:
for xy in XY do
  x,y:=op(xy): z:=x*y/(x + y):
  for a from 2 to 10 do
    for b from 2 to a-1 do
      r:=a/x  + b/y; m:=numer(r); n:=denom(r);
      if n<=10 and m>=2 and m<=n and igcd(a,b,m)=1 
        then N:=N+1; sol[N]:=[x, y, z, a, b, m, n] fi
od od od: time=time[real]()-ti;
'N'=N,entries(sol);

Edit: corrected a<=30 to a<=10. Now 64 solutions in 2 ms.

The radius of convergence can be computed this way:

f:=1/(1-x^2):
x0:=1:
S:=rhs~(`union`(singular(f, x))):
r:=min(seq((abs(x-x0), x in S minus {x0}))):
f=series(f, x=x0); abs(x-x0)<r;

The sphere can be tangent to (some of) the edges in the exterior of the tetrahedron. In this case, Kitonum's conditions must be adjusted.

with(plots):with(plottools):
display(tetrahedron([[0.,0.,0.],[3.,0.,0.],[-4.,3.,0.],[-2.502782406,1.,4.628518959]],transparency=0.2),
        sphere([7.807886553,23.42365966,8.196378356],24.81629404,transparency=0.8,color=yellow),axes=none,
        orientation=[60,70,150]);

not(symbol) in this context is the type of an expression obtained by applying the function not to a symbol.
Examples of this type are not(x), not(Pi).

The type  "not a symbol"  would be Not(symbol).
type("abc", Not(symbol)); # true
type("abc", not(symbol)); # false
 

restart;
H:= k -> < cos(2*k/n*Pi), sin(2*k/n*Pi) >:
n:=100; dt:=0.005; t:=0;
P:=proc(k) global t; Threads:-Sleep(t);t:=t+dt; 
  plots:-display(plot([cos,sin,0..2*Pi]), plots:-arrow(<0,0>,H(k))) end:
Explore(P(k), k=1..n, animate, numframes=n);

Q:=proc(k) 
  plots:-display(plot([cos,sin,0..2*Pi]), plots:-arrow(<0,0>,H(k))) end:
Explore(Q(sqrt(k)), k=1..n^2, animate, numframes=5*n);

This is straightforward (unlike your previous question).

asympt(int((1-t^2)^n, t=0..1),n);

    sqrt(Pi)*sqrt(1/n)/2 + ...

It works for me.

restart;
F:=[-9*a[8]*a[7]*a[9]-5/2*a[3]*a[4]*a[5]*a[9]+6*a[8]^3+5/2*a[3]*a[5]^2*a[8]+3*a[9]
^2*a[6]+3*a[4]^3*a[9]-9*a[4]*a[9]-4*a[4]^2*a[5]*a[8]+3/2*a[4]*a[7]*a[5]^2+9*a[5
]*a[8]-1/2*a[5]^3*a[6], 2*a[3]^4*a[5]-2*a[3]^3*a[4]^2+24*a[3]^2*a[7]^2-48*a[3]*
a[4]*a[6]*a[7]-12*a[3]*a[5]*a[6]^2+36*a[4]^2*a[6]^2+9*a[3]^3-54*a[6]^2, 11*a[5]
^2*a[6]*a[8]+20/3*a[3]*a[5]^2*a[4]^2-12*a[9]*a[4]^2*a[7]-35*a[5]*a[8]^2*a[3]+24
*a[4]^2*a[5]-13/3*a[4]^4*a[5]-90*a[8]^2-11*a[6]*a[4]*a[5]*a[9]+5*a[4]*a[5]*a[7]
*a[8]+10*a[3]*a[7]*a[9]*a[5]+10*a[3]*a[8]*a[9]*a[4]-7/3*a[5]^3*a[3]^2+22*a[8]^2
*a[4]^2+3/2*a[9]^2*a[3]^2-15*a[3]*a[5]^2-3/2*a[7]^2*a[5]^2+36*a[9]*a[7]-27*a[5]
, 10*a[5]*a[3]*a[6]*a[8]+10*a[4]*a[6]*a[7]*a[5]-90*a[7]^2-12*a[6]*a[4]^2*a[8]+
20/3*a[5]*a[3]^2*a[4]^2+11*a[3]^2*a[7]*a[9]-35*a[3]*a[7]^2*a[5]+22*a[7]^2*a[4]^
2+3/2*a[6]^2*a[5]^2-7/3*a[3]^3*a[5]^2-3/2*a[3]^2*a[8]^2+36*a[6]*a[8]+24*a[4]^2*
a[3]-13/3*a[4]^4*a[3]-15*a[5]*a[3]^2+5*a[7]*a[3]*a[4]*a[8]-11*a[3]*a[4]*a[6]*a[
9]-27*a[3], 24*a[3]*a[4]*a[7]*a[9]-87/2*a[3]*a[7]*a[5]*a[8]+24*a[6]*a[4]*a[5]*a
[8]+33/2*a[3]*a[6]*a[5]*a[9]+54*a[4]^3-9*a[4]^5+9/2*a[8]*a[3]^2*a[9]+51*a[4]^2*
a[7]*a[8]+9/2*a[7]*a[5]^2*a[6]-7*a[3]^2*a[5]^2*a[4]+16*a[4]^3*a[5]*a[3]-45*a[6]
*a[4]^2*a[9]-18*a[3]*a[4]*a[8]^2-45*a[5]*a[3]*a[4]-18*a[4]*a[7]^2*a[5]-81*a[4]+
81*a[6]*a[9]-135*a[7]*a[8], 138*a[6]*a[7]*a[4]*a[9]-129*a[3]*a[4]^2*a[7]*a[5]-\
147*a[6]*a[7]*a[5]*a[8]+187/2*a[3]^2*a[4]*a[5]*a[8]-360*a[4]^2*a[7]+62*a[7]*a[4
]^4+27*a[3]^2*a[9]+18*a[7]^3*a[5]-48*a[3]*a[6]*a[8]*a[9]-129/2*a[3]*a[6]*a[5]^2
*a[4]-243*a[6]*a[4]*a[5]+85*a[4]^3*a[5]*a[6]-21*a[3]^2*a[4]^2*a[9]+288*a[8]*a[3
]*a[4]+288*a[5]*a[3]*a[7]+91/2*a[7]*a[5]^2*a[3]^2-76*a[4]^3*a[8]*a[3]-24*a[3]*a
[7]^2*a[9]+192*a[3]*a[7]*a[8]^2-132*a[6]*a[4]*a[8]^2+33*a[6]^2*a[5]*a[9]+9/2*a[
3]^3*a[5]*a[9]-30*a[7]^2*a[4]*a[8]+486*a[7], 270*a[6]^2*a[4]*a[9]-3*a[3]^2*a[4]
^2*a[8]-33/2*a[3]^3*a[4]*a[9]+57/2*a[3]^3*a[5]*a[8]+261*a[3]*a[7]^2*a[8]-72*a[6
]^2*a[5]*a[8]+621*a[7]*a[3]*a[4]+54*a[3]*a[6]*a[8]^2+405*a[5]*a[3]*a[6]+87/2*a[
3]^2*a[5]^2*a[6]-219*a[4]^3*a[7]*a[3]-1134*a[6]*a[4]^2+270*a[6]*a[4]^4+108*a[3]
^2*a[8]+108*a[4]*a[7]^3+363/2*a[4]*a[3]^2*a[7]*a[5]-171*a[3]*a[6]*a[7]*a[9]-450
*a[6]*a[4]*a[7]*a[8]-285*a[6]*a[4]^2*a[5]*a[3]+972*a[6]]:
A:=indets(F)[];
X:=seq(x||i,i=3..9);
G:=Groebner:-Basis(F,tdeg(A));
FX:=eval(F, [A=~X]):
GX:=Groebner:-Basis(FX,tdeg(X));
evalb( eval(G, [A=~X]) = GX);  # true

BTW, I don't recommend using the document mode (and 2D input) for programming. For example it is hard to run your worksheet line by line.

As mmcdara said, this can be computed using the Laplace method. Unfortunaly he made a mistake, the maximum point found was outside the interval of integration.

restart;
J:=  Int(exp(-n*(x*cosh(t)+t)), t = 0 .. infinity):
h:= convert(series(-n*(x*cosh(t)+t),t,3), polynom):
int(exp(h),t=0..infinity) assuming x>0,n>0:
A:=asympt(%,n,2) assuming x>0;

     A := (1/n+O(1/n^2))/exp(n*x)

evalf( eval([J=A], [x=1, n=100]));  # check
evalf( eval([J=A], [x=1/4, n=100]));

     [3.683935796*10^(-46) = 3.720075976*10^(-46)+3.720075976*10^(-44)*O(1/10000)]
     [1.385347786*10^(-13) = 1.388794386*10^(-13)+1.388794386*10^(-11)*O(1/10000)]

For simplicity, I will suppose that f : (0, oo) --> R  is continuous and the relation

Int(f(x), x=a..b) = Int(f(x), x=1/b..1/a)    (*)

holds for all 0 < a < b < oo.

(*) is equivalent to 

Int(f(x) - f(1/x)/x^2 , x=a..b) = 0, for all a,b>0
so, f(x) - f(1/x)/x^2 = 0, for all x>0.

The general form for f is given by:

f(x) = piecewise(x<1, u(x), u(1/x)/x^2),

where u : (0,1] --> R is an arbitrary continuous function.

Note. Your example f(x) = 1/(x^2+1) satisfies f(1/x^2)/x^2 = f(x).