Rewrite the recurrence:

rsolve({f(n+1)-f(n) = 1+f(n), f(1)=1}, f(n))

      2^n-1

Not all mathematical functions have evalhf implementation. pochhammer is one of them. See ?evalhf/fcnlist

You will have to use evalf instead.

F := (x,y,z) -> sin(1/x+1/(y*z^2)):
eval(series(F(x/t,y/t,z/t),t), t=1);

     1/x + 1/(y*z^2) - 1/(6*x^3) - 1/(2*y*z^2*x^2) + 1/(120*x^5)

You can use the package LinearAlgebra:-Generic.
It has not a Rank command, but there is ReducedRowEchelonForm which can be used easily to compute the rank.

If you want to check the inequality, use a random test.
Probably it's easier to construct your own:

restart;
fn := (a^(n+1)-b^(n+1))< c*(a^n-b^n);
r:=rand(0.0 .. 1):
rn:=rand(1..10):
N:=100000: OK:=0:  ex:=[];
to N do
  c:=r(): b:=c*r(): a:=b*r(); n:=rn():
  if (fn) then OK:=OK+1 else ex:=[a,b,c,n] fi
od:
'OK'=OK,'N'=N, 'counterexample'=ex;

    OK = 7393, N = 100000, counterexample = [0.2282758836, 0.6582892290, 0.6718939194, 5]

So, it is not true, and you have a counter-example.
(It is possible to exit the loop as soon as you find one.)

Download linsys.mw

Just delete the # in the last line and then press ENTER.

simplify(expand(f)) assuming x>0;

So, you want an example. The problem is interesting.

For any a[n]>0 converging to 0 very fast, P(n) will have n (negative) roots.
A concrete example is a[n] = 1 / 2^(3^n).

Check:

P:= n -> add(x^k/ 2^(3^k), k=0..n):
seq( [n, nops([fsolve(P(n))])], n=1..14);

    [1, 1], [2, 2], [3, 3], [4, 4], [5, 5], [6, 6], [7, 7], [8, 8], [9, 9], [10, 10], [11, 11], [12, 12], [13, 13], [14, 14]

Unfortunately for n>14 Maple will have problems here because 2^(3^15) has > 4*10^6 digits and the roots are also huge.

The integral can be computed.

J := fun11[1]
result11:=eval(subsop(2=s,J));

You have now the indefinite integral. It is a huge elementary expression (because you have a lot of parameters).
I do not know why you need it symbolically.
Anyway, in order to obtain the definite integral it's not possible to use directly the Newton-Leibniz formula.
(You can of course do it but the result will be generic and twice bigger).

You will have to check that  result11 is continuous in your interval; this will depend on the parameters.
In case of discontinuities, the integral must be split first. That's why Maple does not do it: there are too many possibilities, and the expressions involved are huge.

 

restart;
with(LinearAlgebra):
A := Matrix(4, 4, {(1, 1) = 1, (1, 3) = -1, (1, 4) = 3, (2, 2) = 2, (2, 3) = 1, (3, 1) = -1, (3, 2) = 1, (3, 3) = 6, (3, 4) = -1, (4, 1) = 3, (4, 3) = -1, (4, 4) = 10}, fill = 0):
b := Matrix(4, 1, {(1, 1) = 0, (2, 1) = -2, (3, 1) = -1, (4, 1) = -1}):
x := Matrix([[x1], [x2], [x3], [x4]]):
f := x -> (((1/2*Transpose(x)) . A) . x) + ((Transpose(b)) . x):
v0 := Matrix([[0.], [1.], [0.], [0.]]):
g := x -> (A . x) + b:
while   (Norm(g(v0)) > 1e-6) do;
  alpha0 := solve(diff(f(v0 - g(v0)*alpha)[1, 1], alpha) = 0);
  v0 := v0 - alpha0*g(v0);
od:
v0, g(v0);

Note that LinearSolve(A,b);  gives

f:=t -> t^2;
a := n -> 1/Pi*int(cos(n*t)*f(t),t=0..2*Pi):
b := n -> 1/Pi*int(sin(n*t)*f(t),t=0..2*Pi):
S := (n,t) -> a(0)/2+add(a(k)*cos(k*t)+b(k)*sin(k*t),k=1..n):
plot([f(t),S(8,t)], t=0..2*Pi);

It can be shown easily using the initial form that the function is not periodic (so, the "period" is infinity).
It is however almost periodic (see wiki) .

For x > 0:

F1(x) = -c*x^2 - x*sin(x^2) - 2*x + 1,  F2(x) = 2/sqrt(x) + c;

For x<0, F2(x) is arbitrary.

(I suppose that Phi = 1 only for x=y, the problem being unclear wrt this).

It works for me in Maple 2018, 64 bit.

For n=9 you already have 261080 graphs!
In combinatorics you should always estimate the size of the problem and don't ask for the impossible.