@acer 

Just want to mention that in this case allvalues does not garantee that all the critical points are returned,
as _ValuesMayBeLost shows.

@Carl Love 

I prefer the new Iterator

f := proc(u)
     local M:=convert(u,list);
     `if`(LinearAlgebra:-Determinant(Matrix(4,M)) mod 2 <>0, cat(M[]), NULL)
     end;
P:=Iterator:-CartesianProduct([seq(0..1)]$16):
GROUP:=seq(f(M),M=P):   
nops([GROUP]);

It is fast enough and straightforward. In other cases (when conversions are not needed) its speed is much better compared with other methods.
(Anyway the speed here is irrelevant since bigger dimensions are out of the question.)

@Carl Love 

LinearAlgebra:-Determinant works as well and is twice faster on my computer.
(Det would be better only for a big PRIME but in this case the enumeration is impossible.)

@one man 

With your prefered method (where x,y are supposed to be Re and Im for z) ==> complex values for x,y.
Are you really satisfied with this? Note that Maple uses complex numbers by default.

@acer 

Is seems that you found a mistake in the code, but this would solve only the first bug (and partially the third, because this way at least `infinity` is replaced by an unevaluated `minimize`).

But Maple should be able to compute all three examples, because it can find the critical points:
f:=(x+y)^2+cos(y):
s:=solve([diff(f,x),diff(f,y)],[x,y], allsolutions);
      s := [[x = -Pi*_Z1~, y = Pi*_Z1~]]
eval(f,s[]);
      (-1)^_Z1~

The OP probably does not realize that  w=f(r), z=g(r)  is a parametric representation of a curve in the (w,z) space and its graph it's what he wants.

Example:   w=r+sin(6*r), z=r+sin(7*r),   r in [0,2*Pi];

plot([r+sin(6*r), r+sin(7*r) ,r=0..2*Pi], labels=["w","z"]);

@Markiyan Hirnyk 

Because the polynomials are very simple. But try this for

n^17-n^12-n^5+7
and
a^17-85*a^16+3400*a^15-85000*a^14+1487500*a^13-19337501*a^12+193375060*a^11-1519376650*a^10+9496121250*a^9-47480778125*a^8+189924350000*a^7-604311312500*a^6+1510804062499*a^5-2905466796850*a^4+4150820312250*a^3-4151035155000*a^2+2594580075000*a-763183590618

Here it is easier by hand!

@Markiyan Hirnyk 

If in Mathematica one wants to retrive Kitonum's solutions, a similar by hand substitution will be needed.

@Markiyan Hirnyk 

You are right. I have corrected the second example, which was answered correctly by Maple.

@Markiyan Hirnyk 

Workarounds? But these are bugs too!

@tomleslie 

Yes, it is Pi^2/8 with "probability 1".

It reduces to
J:=2*Int(arcsinh(cosh(x))-x, x=0..infinity);
identify(evalf(J));  #Pi^2/8

is(value(J)=Pi^2/8); #false

It seems to be hard to compute J symbolically.

@John Fredsted 

Theoretically it could be formulated as sequence of optimization problems:
find the vectors vi = [xi1,xi2,...,xi8] (i=1,...,n) satisfying the constraints ||vi||=1, ||vi-vj|| >= 1.

Then increase n until a solution does not exist.

But practically the number of the constraints n*(n+1)/2 is probably too large for n>50 say, to hope for a solution.

@John Fredsted 

Why do you say that? A double-check is always better than none!

Edit. I mean e.g. to use Maple to find a maximal set of vectors V on the sphere ||x||=r in R^8 such that ||x-y|| >= r, for x,y in V, x<>y.

@John Fredsted 

A simple inspection is enough. E.g. if both x,y have integer entries then x-y can have 1,2,3 or 4 nonzero entries, and in the case of a single nonzero entry, this must be +/-2. So ||x-y|| >= sqrt(2) [attained if x-y has 2 nonzero entries +/-1].

P.S. It would be nice to use Maple to find such a configuration.

I don't know much about sphere packing, but I don't understand why do you need a Maple check here, when mathematically this is almost obvious: for E8, if x,y are distinct "root" vectors, ||x||=sqrt(2), ||x-y|| >= sqrt(2).
Probably in one of the Leech constructions the situation is similar.