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

{seq(msolve({a = k, a*b = 2391}, 10000), k = 1 .. 9999)};

with(ImageTools): 
img0:=Read("smile.jpg"):   str:="Just a test!":
img1:=Scale(img0,1..50):
img:=Create(400,400,channels=3,background=white):
simg:=(x,y) -> SetSubImage(img,img1, min(max(floor(400-400*y-25),1),350), min(max(floor(400*x-25),1),350)):
Explore( plot(x^2,x=0..1, background=simg(a,a^2), axes=none, title=str[1..floor(15*a)]), a=0..1., animate,loop,autorun);

(The file "smile.jpg"  must be in the current directory).

Your system has a sigularity at t=0.
Just replace in ICS:  S(0) = 30 with say S(0.001) = 30  etc. It will work.
 

It was answered recently here.

You cannot use the operator * as a name. So, use  beta[`*`], or even better `beta__*`.

Also, remove the last comma in params.
 

In Maple there are many options for numerical integration (including compiled external libraries) which usually can increase considerably the speed. You should choose the proper ones for your case. Read the help pages:

?evalf,int

and also the provided links.