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

The shape of a quadrilateral inscribed in the unit circle and having an inscribed circle.

f := (a,b) -> 1/2*(b+arccos(sin(2*a-b)/tan(b))):
dis:=(A,B)-> sqrt(inner(A-B,A-B)):
bisA:=proc(A,B,C)
  local b,c,M;
  b,c:=dis(A,C),dis(A,B);
  M:=(b*B+c*C)/(b+c);
  x*(A[2]-M[2])+y*(M[1]-A[1]) + A[1]*M[2]-A[2]*M[1]
end:
P:=proc(a0,b0)
  local a:=evalf(a0), b:=evalf(b0), c:=f(a,b), p1,p2,p3,p4, r, II;
  if a0>b0-0.001 then return plot() fi;
  p1,p2,p3,p4:=[1,0],[cos(4*a),sin(4*a)],[cos(4*b),sin(4*b)],[cos(4*c),sin(4*c)];
  II:= eval([x,y], solve( [bisA(p2,p1,p3), bisA(p3,p2,p4)]) );
  r:=evalf(eval(dis(II,p2)*cos(b))); 
  plots:-display(plot([p1,p2,p3,p4,p1]), plottools:-circle(II,r, color=blue), plottools:-circle([0,0],1));
end:
Explore(P(a,b), parameters=[a=0.01..Pi/2-0.02, b=0.02..Pi/2-0.01], initialvalues=[a=Pi/8,b=Pi/4]);

###    0 < a < b < Pi/2