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.

f must be monotonic and unbounded.
Computing diff(f,x) it follows that   k∈(-∞,0)∪[3/2,+∞)

No need for a "real" algorithm; a double loop is enough.
(unless you are looking for best speed, but the gain would be small).

restart;
# 1 < h < k <= n
n:=30:
a:=1/1000:
X:=sort(LinearAlgebra:-RandomVector(n, generator=rand(0..10^3)));

t:= (h,k) -> min(X[k]-X[k-1], (X[h+1]-X[h])/(n+1-k)) * (n+1-k)*(n+k-2*h):
A := Matrix(n, (h,k)->`if`(h<k and h>1, t(h,k), -infinity)):
h,k:=max[index](A):
`if`( a*n*add(X) <= A[h,k], ['h'=h,'k'=k], FAIL);

The first two (iterated) limits do not exist, because the inner limits do not exist.
Maple answers correctly:
f:=(x+y)*sin(1/x)*sin(1/y);
limit( f, x=0); # ==> an interval

The last limit is 0 (provided that x,y are real; the limit does not exist in C);

limit(f, {x=0, y=0});  # ==> 0

Yes, it's a bug.

A := <"11",12,13;21,"22",23>:
DocumentTools[Tabulate](A, width=40, fillcolor=((x,i,j)->`if`(irem(j,2)=1,cyan,red)) ):  #ok

DocumentTools[Tabulate](A, width=40,      color=((x,i,j)->`if`(irem(j,2)=1,cyan,red)) ):  # only the strings are colored

Here is the direct approach, i.e. defining the procedures.

P := proc(t) option remember; a*ED(t - 1) + P(t - 1)   end;
ED := proc(t) option remember; DC(t) + DF(t) end;
DC := proc(t) option remember; c*(P(t) - P(t - 1)) end;
DF := proc(t) option remember; b*(F - P(t)) end;
a := 1;
c := 0.75;
b := 0.2;
F := 100;
P(0) := F;
P(1) := F + 1;

P(10), DF(10);

     99.38828842, 0.122342316

plot( [seq([t,P(t)],t=0..100)]);  # you may want to add style=point

In principle, Student:-Precalculus:-CompleteSquare should work but in this case it complicates things.
I prefer to use my generalized version:

SQR:=proc(P::polynom(anything,x), x::name)
local n:=degree(P,x)/2, q,r,Q,R,k;
if not(type(n,posint)) then error "degree(P) must be even" fi;
Q:=add(q[k]*x^k,k=0..n-1) + sqrt(lcoeff(P,x))*x^n;
R:=add(r[k]*x^k,k=0..n-1); 
solve({coeffs(expand(P-Q^2-R),x)}, {seq(q[k],k=0..n-1),seq(r[k],k=0..n-1)});
eval(Q,%)^2+eval(R,%)
end:

expr1:=-lambda-(1/2)*kappa__c-gamma__p-(1/2)*sqrt(-16*N*g^2+4*lambda^2-8*lambda*gamma__p+4*lambda*kappa__c+4*gamma__p^2-4*gamma__p*kappa__c+kappa__c^2):

evalindets(expr1, sqrt, t -> SQR( op(1,t), lambda)^op(2,t));