Of course.
The mapping (x,y) -> `if`(x+y<1, [x,y],[1-x,1-y])   applies the square [0,1]^2 onto a standard triangle similar to the mentioned Dirichlet distribution.

@Markiyan Hirnyk 

My code does not use incidence verification. It constructs effectively area preserving mappings just like Mathematica.

@Markiyan Hirnyk 

@Heitor 

It was just a joke, probably not a good one.

@torabi 

A general advise (at least for the future): it is a good idea to present mathematically the problem.  Otherwise what alternative method are you hoping for? It would imply a "reverse engineering" of your computations and probably nobody will want to spend his time for this. E.g. probably Runge-Kutta makes sense for the original problem, but not for the system you have presented.

@ecterrab 

What is the performance impact of Assume? After Physics:-Assume(...), Maple seems to use more memory.
It would be interesting to know why this operating mode was not included in the standard assume command (maybe via a kernelopts(assume) switch).

If you want to suppress the output for .txt use colon (:) after read
see also ?interface,echo

To use assumptions in .m file save `property/object` as mentioned in the answer.

@Kitonum 

Please try the version:

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, compile = false):
GROUP:=seq(f(M),M=P):   
nops([GROUP]);

which disables the compilation, and see if it works for you. It is just a little slower.

But you should try to solve the compiler problem, see ?compiler

@Carl Love 

GL(2,64) has > 1.6*10^7 elements.

But note that the methods here produce incorrect results for a non-prime q.

@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!