You have only two variables (eta and u). So, what kind of 3D do you mean?

The legends are separated from the plot area in Maple. You may use textplot etc if you really want this.

@mmcdara Yes, I also use this often. Vote up.
Unfortunately it may fail even for simple expressions.

J:=Int(Heaviside(1-x^2-y^2)*x^2, [x=-1..1, y=-1..1]):
K:=Int(Heaviside(1-x^2-y^2)*x^2, [x=-infinity..infinity, y=-infinity..infinity]):
value(J) = evalf(J);
#                        0 = 0.7853981634
identify(%);
#                        0 = Pi/4
value(K);
#                          undefined
evalf(K);
#                          0.7853981634

Note that even int(Heaviside(1-x^2-y^2), [x=-1..1, y=-1..1])  fails (it gives 0 too).

@mmcdara It was approx 70 sec,  almost all for the second eliminate (for I2).

@mmcdara  If `union`(S) = {a1,...,an}  and  nops(S) >> n,  probably a good candidate for F(S) is  { {}, {a1}, ..., {an} }.

@tomleslie I don't see other solution if the present functionality is maintained. The default labels could be `1`, `2`,... or maybe 1., 2., ... (non-integer). Or, it may remain as it is, with a warning.

@Earl If the contact points are allowed to be anywere on the lines, any conic is possible.

@Carl Love I ment f (z,A,phi) has values arbitraryly close to 0 in any nbd of complex infinity (complex infinity being an essential singularity of f(., A, phi)). So, it will be impossible to find numerically all the roots, if a bounded region is not known.

You will need a theoretical analysis first because:
1. solve may omit solutions
2. f (z,A,phi) has values arbitraryly close to 0 (by Casorati–Weierstrass theorem).

BTW. It's not a good idea to use Zeta as variable (it's reserved for Riemann's function). 

@Carl Love  It would be interesting a 10 minutes contest: given Count, explain what it does and how it works.

.maple files are not accepted by mapleprimes. Use .mw or compress into a .zip.

@Rouben Rostamian  Yes, it's a vector graphics, but I suspect that it was produced from a bitmap.
In Maxima the EPS is obtained in a second and has 200KB.

@Rouben Rostamian  In Windows, exporting 3d graphics to EPS does not work well at all.
For example 

p:=plot3d([x^2+y^2,1-x^2-y^2], x=-1..1,y=-1..1);
plottools:-exportplot("d:/tmp/aatest.eps", p);

 needs several minutes, has 12MB and artifacts.
It is probably a vectorized version of a bitmap. Useless.

@mmcdara In Maple 2020+

sort(1+x,descending):
# ...
g := x -> 1/(1+x):
(g@@3)(x);

@Carl Love You wrote  "Because the module references are linked into the code at the time that the code is read, not the time that it's executed".

But when uses IM is parsed, IM was already parsed, so the exports are known and it would be enough for the parser to "observe" that add_three is exported by IM and to add the IM:- prefix.
I don't see any reason to reject uses IM.

@mmcdara 

1. I think that this can be set with the Typesetting package, e.g. with the interactive Typesetting[RuleAssistant]()
I do not use it much. Anyway, after with(Logic) the desired typesetting is automatic.

2. A name beginning with & is treated as a binary operator (or prefix unary). See ?neutral operators