My impression is that MultiSeries:-limit was not designed for complex, and this option should be removed (just use the usual :-limit(...,complex) which seems to work rather well).
All next answers are strange:

with(MultiSeries):
limit(sin(1/z)*z,z=0, complex);
                               0
limit(sin(I/z)*z,z=0, complex); #ok but lucky
                           undefined
limit(GAMMA(z),z=0, complex);
                            infinity
limit(1/z, z=0, complex);
                            infinity
limit(-1/z, z=0, complex);
                           -infinity
limit(I/z, z=0, complex);
                           infinity I
limit(-2*I/z,z=0, complex);
                          -infinity I

The Lebesgue-Schoenberg version works easily in higher dimensions:
(For Hilbert/Peano extra work is needed.)

F:=proc(t)
   local u:=abs(2*round(t/2)-t);
   piecewise(u<=1/3,0, u<2/3, 3*u-1,1)
end:
PN:=4:
Pd:=3:  # number of dimensions
for i to Pd do
 P[i] := unapply( add(F(3^(Pd*k+i-1)*t)/2^(k+1),k=0..PN), t):
od:

plots:-spacecurve([P[1](t),P[2](t),P[3](t)],t=0..1, numpoints=100000);

@Carl Love 

But then e.g. a symbolic integration is out of the question.

P.S. I agree that this is not a bug sensu stricto but Maple should do much better in such situations.

@rlopez 

Isn't 'singsol' = all a more natural option?

@Preben Alsholm 

The function is, is also used (outside the assume facility) as a "smarter" boolean evaluator.
E.g. evalb( 10^(1/3)<3 ); versus is( 10^(1/3)<3 );
I don't see why including elif in the definition would be useful, except that g(x) is now undefined for a symbolic x.
I can understand why Maple gives these results (due to evaluation rules) but I still think that there exist much better solutions (for int and plot e.g) in such situations.

@mostafajani 

This happens due to the evaluation rules in Maple.
Note that  f(x)  returns x^2  because  is(x<0)  returns false (when x is symbolic).
Now, when plot evaluates its argument it gets x^2 and plots it.

Use:
plot('f'(x), x = -2 .. 2);
or
plot(f, -2 .. 2);

@tomleslie 

p := x^4-17001090662707204*x^2+4:
s1:=select(u -> is(u>0), [solve(p,explicit)]);
    s1 := [46099201*sqrt(2)-37639840*sqrt(3), 46099201*sqrt(2)+37639840*sqrt(3)]
evalf(s1); evalf[30](s1);
    [-0.3e-1, 1.303882305*10^8]
    [1.53388077417338*10^(-8), 1.30388230537526675954819973540*10^8]

# So, s1 is correct; the first solution is >0 but evalf needs >10 digits to see this

s2:=select(type, allvalues~([solve(p)]), positive &under evalf);
    s2 := [-46099201*sqrt(2)+37639840*sqrt(3), 46099201*sqrt(2)+37639840*sqrt(3)]
evalf(s2); evalf[30](s2);
[0.3e-1, 1.303882305*10^8]
[-1.53388077417338*10^(-8), 1.30388230537526675954819973540*10^8]
#s2 is not correct

s3:=solve({p, x>0});
    s3 := {x = RootOf(_Z^4-17001090662707204*_Z^2+4, 0. .. 0.625e-1)}, {x = 1.303882305*10^8}
evalf(s3); evalf[30](s3);
    {x = 1.533880774*10^(-8)}, {x = 1.303882305*10^8}
    {x = 1.53388077417339100429211486174*10^(-8)}, {x = 1.303882305*10^8}

# So, s3 is correct
# It seems that I was wrong and only Carl's solutions contains negative numbers.
# This is not a bug though because s2 selects the positive solutions using the current Digits (=10)
# while s1 and s3 are obtained internally and Digits is most probably increased.
# [actually I was surprised in the first post when it seemed that s3 <> s1]
#
#  I hope that your head is OK now  :-)

Actually both slick/non slick methods fail for the polynomial
p := x^4-17001090662707204*x^2+4;
I'd suggest:

select(u -> is(u>0), [solve(p,explicit)]);

which seems to work without modifying Digits.

@Zeineb 

I said that |f(z)| < 1 at the corners i.e. for x,y in {- (n+1/2)*Pi , (n+1/2)*Pi },  n in Z. And this is true.

The general situation is:

|f(z)| < 1  if and only if  

x in R  and  n*Pi + Pi/4 < y < n*Pi + 3*Pi/4  for some integer n.

@Kitonum 
The OP wanted this way. Otherwise is even simpler:
GetCoeff:=proc(f,u, x)  # x = the set of vars
  if u=0 or u=1 then return eval(f,x=~0) fi;
  eval(foldl(coeff,f,op(u)), x=~0);
end proc:

@Rouben Rostamian  

eval( m, [x=Pi*(n+1/2),y=Pi*(n+1/2)]):  # it's enough, cos,cosh being even
simplify(%) assuming n::integer:
convert(%,tanh);
              tanh((1/2)*Pi*(2*n+1))^2
So, <1.

@Kitonum 

There is another problem. If C is the cuboid [a1,a2] x [b1,b2] x [c1,c2], this approach plots f( boundary(C) ) , and this is not necessarily the boundary of f(C),  which is supposed to be plotted.

@Markiyan Hirnyk 

int( sin(x^2 - y^2), dxdy )  = 0

(by symmetry, the integrand is odd)

@toandhsp 

You may produce a warning (or an error if you prefer). E.g.

eq:=proc(L::listlist)
  local e:=primpart(LinearAlgebra:-Determinant(<Matrix([[x,y,z],L[]])|<1,1,1,1>>));
  if e=0 then WARNING("Collinear points: %1",L) end if;
  e*signum(lcoeff(e))=0
end proc:

@toandhsp 

eq:=proc(L::listlist,x::name,y::name,z::name)
  primpart(LinearAlgebra:-Determinant(<Matrix([[x,y,z],L[]])|<1,1,1,1>>));
  %*signum(lcoeff(%))=0
end:

Note that this also works if the list contains floats.