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

{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.