There are essentially two methods.

1. Use DirectSearch.
GlobalOptima(f(x,y), {x=3..4, y=0..1}, maximize);
The approx max is (almost) 0.

2. Use calculus.
Solve the system [diff(f(x,y),x), diff(f(x,y),y)].
Maple finds a unique solution in the domain: x=4, y=0, ( f(4,0)=0 ).
It remains to compute the max on the boundary: f(x,0), f(x,1), f(3,y), f(4,y)
(it's easy).

Maple has not a package for solving general functional equations.
A few years ago I saw at a conference a group of hugarian mathematicians developping such a package, but I cannot find the references (they should be somewhere on my computer).
Anyway, the equations you mention do not seem to be difficult. The general solution for the first one (for x>0) is:

F(x) = piecewise( x<1, a(x), x^(-r)*a(1/x) )

where a(x) is an arbitrary function for 0 < x <= 1.

@escorpsy 

Your integral diverges. It has a singularity at Zeta=Pi/2. You may check with MultiSeries:-series(phi*W, Zeta = Pi/2)

or compute some values around Pi/2.

P.S. You should not use Zeta as a variable, it is the well known Riemann's function; use zeta instead!

The irreducibility is a very difficult math problem and needs a precise context.
Let's consider the case of polynomials in Z[x,y], which are linear in parameter d (integer too).

restart; # example
randomize():
expand((7+randpoly([x,y]))*(66+randpoly([x,y]))): eval(%,coeffs(%)[5]=d):
p:=%;

p := d*x^7*y^2 - 374*x^9*y - 2510*x^7*y^3 - 1224*x^5*y^5 - 1156*x^9 - 2176*x^5*y^4 + 825*x^4*y^5 + 5100*x^2*y^7 - 2856*x^8 - 1088*x^7*y - 1512*x^6*y^2 + 2550*x^4*y^4 - 3249*x^6*y - 3820*x^4*y^3 + 6300*x^3*y^4 - 2176*x^2*y^5 - 476*x^6 - 544*x^5*y - 422*x^4*y^2 - 288*x^3*y^3 - 4522*x^5 + 1910*x^4*y - 3326*x^3*y^2 + 4488*x^2*y^3 - 2775*x*y^4 + 1200*y^5 + 2244*x^4 - 224*x^3*y + 525*y^4 + 7431*x^3 - 816*x^2*y + 1184*x*y^2 - 512*y^3 - 357*x^2 - 224*y^2 - 2442*x + 1056*y + 462

Here we can try brute force

n:=10000:
for d in {seq(-n..n)} do
  if not irreduc(p) then print('d'=d); break fi od:

                        d = -964

factor(p);

      -(11*x^4*y + 68*x^2*y^3 + 34*x^4 + 84*x^3 - 37*x + 16*y + 7)*(34*x^5 + 18*x^3*y^2 + 32*x^3*y - 75*y^4 + 51*x^2 + 32*y^2 - 66)

You should be pleased about how  '*' and '/' work for inequalities, because this way you will avoid incorrect results.
E.g. try to multiply   -2 < -1   and   1 < 3.

restart;
with(GraphTheory):
with(SpecialGraphs):
K4 := CompleteGraph(4):
IsPlanar(K4, 'F'):   # true

G:=K4:
n:=max(Vertices(G)):
V1:=[Vertices(G)[],n+1]:
E1:=Edges(G) union {seq({i,n+1},i=1..n)}:
G1:=Graph(V1,E1):
IsPlanar(G1, 'F1'); # iff G is Outerplanar

        false

I had a similar problem with SMTLIB:-Satisfy in Maple 2018 (on Windows 7).

Installing the Visual C++ Redistributable for Visual Studio 2012 Update 4

https://www.microsoft.com/en-us/download/details.aspx?id=30679#

solved the problem.

Maple does not miss solutions here. But solve uses the (documented) convention that the arguments of abs are assumed to be real.

Note that the equation has infinitely many complex solutions; they can be found in terms of RootOfs (or even explicitely, but with huge expressions).

However, in this case Maple assumes that the unknown itself is real, not only the argument of abs. (There are 4 complex solutions having the argument of abs equal to -1).

It's the Java GUI. It happens sometimes. I miss the Classic interface. Java could be nice for developers, but for users ... slow and buggy.

Here is a general solution.
Probably it cannot be accepted as a homework solution, but it's a good exercise for you to understand how it works (Maple+maths).

n:=4: # take n:=5 for your problem
f := ( 1 + a * mul(x-k, k=1..n) )/x:
f := unapply(simplify(eval(f,solve(eval(numer(f),x=0),{a}))), x):
seq( 'f'(k)=f(k), k=-1..7 );

       f(-1) = 4, f(0) = 25/12, f(1) = 1, f(2) = 1/2, f(3) = 1/3, f(4) = 1/4, f(5) = 0, f(6) = -2/3, f(7) = -2

For n = 34, R[n] alone is a rational number whose denominator has about 5*10^9 digits!

Are you sure you really want this?

You can use the explicit option in solve and avoid allvalues.

F := [x^2+y+z-1, y^2+x+z-1, z^2+x+y-1]:
Sol:=solve(F, explicit);
    Sol := {x = 0, y = 0, z = 1}, {x = 0, y = 1, z = 0}, {x = 1, y = 0, z = 0}, {x = 2^(1/2)-1, y = 2^(1/2)-1, z = 2^(1/2)-1}, {x = -1-2^(1/2), y = -1-2^(1/2), z = -1-2^(1/2)}
seq(simplify(eval(F,s)), s=Sol);
     [0, 0, 0], [0, 0, 0], [0, 0, 0], [0, 0, 0], [0, 0, 0]
# or ...
seq(simplify(subs(s,F)), s=Sol);
    [0, 0, 0], [0, 0, 0], [0, 0, 0], [0, 0, 0], [0, 0, 0]

If we insist on using LSystem, a more natural solution is:

restart;
with(Fractals:-LSystem): with(plots):
state:= "F-F-F":
cons:= proc(n)
  local x:=sprintf("%.10f",((1/3.)^n));
  [ "F"="draw:"||x, "+"="turn:-60", "-"="turn:-240" ]
end:
rules:= ["F"="F+F-F+F"]:
display(seq(LSystemPlot(Iterate(state, rules, n), cons(n))$4, n=0..6), insequence);

One obtains the same animation as in my first direct solution.

nextkoch:=(a,b)->
((2*a+b)/3, (2*a+b)/3+ (b-a)/3*evalf(exp(-Pi/3*I)), (a+2*b)/3):
L:=evalf([0,1,1/2+I*sqrt(3)/2,0]): N:=6: LL[0]:=L:
for j to N do
  L:=[ seq( [ L[i], nextkoch(L[i],L[i+1]) ][], i=1..nops(L)-1), L[-1] ]:
  LL[j]:=L;
od:
plots:-animate( plots:-complexplot, [LL[floor(n)]], n=0..N+.99, 
  scaling=constrained, axes=none, frames=4*(N+1), paraminfo=false);