@arjangash If you really have no problem about optimization, it should be obvious.
Here is a simpler example:

restart;
Z:=(a,b,g)->a+b+g:
constr:=(a,b,g) -> (a>=0,b>=0, a^2+b^2-2*a+1<=g):
Zmax:= g -> Optimization:-Maximize(Z(a,b,g), [constr(a,b,g)])[1]:
plot(Zmax, 0.1 .. 2.2);

@acer 

You are right, I simply forgot about having reported it, sorry! Should I delete the question?

In my opinion the problem is not to find workarounds (which are not complicated) but to understand where exactly is Maple's fault.

Maple can be used to compute the limit of the sequence

n -> sqrt(n) * (sin@@n)(1);

Can MMA compute it?

@Markiyan Hirnyk 

It would have been nice if MMA had given the simplified result sqrt(2)/4.
Probably it simply uses the fact that the sequence exp(I*n) is dense in the unit circle.

Maybe some day we will get something like:
f := sin(n)/(3 + cos(n)):
limsup[discrete](f,n)
     sqrt(2)/4;

For the moment:
maximize(f, n=0..2*Pi);
    
 

Just a short remark.
If the sequence is given by a procedure f, then its limit is L iff
limit( f(floor(x)), x=infinity) = L.

This way, the assumption is not necessary.

Unfortunately this does not work because floor has a very poor integration in Maple [I don't know why].
A simple (and embarrassing) example is:

simplify( floor(x) ) assuming x>3,x<4;
    floor(x)

@Markiyan Hirnyk 

It is OK for your surface defined by x^2+y^2+z^2=1, x^2-z^2<=1
but not for Mariusz' which was the boundary of the body  x^2+y^2+z^2<=1, x^2-z^2<=1.
In the previous comment both situations were answered.
But indeed in Mathematica code I see an "<"; probably this means the exclusion of this part of the boundary [I do not know the MMA syntax for this].

@Mariusz Iwaniuk

Download surface-int.mw

@digerdiga 

For polynomials you should use:

restart;
q:=(x-1)*(x-2)*(x-3):
sol:=c -> fsolve(q+c,x,complex):
r:=c -> max(Re~(select(t->abs(Im(t))<1e-10, [sol(c)]))):
plot(r, -8..8, discont, color=green);

Or, better:

restart;
q:=(x-1)*(x-2)*(x-3):
sol:=c -> fsolve(q+c,x):
r:=c -> max(sol(c)):
plot(r, -8..8, discont, color=green);

@Kitonum 

I am sorry that you, the proposer of the test, have trouble in getting the answer.
It seems that SolveTools:-SemiAlgebraic has a bug in the 32 bit version [or, it could be a memory management problem].
On a 32 bit computer I also get the same error for:

SolveTools:-SemiAlgebraic(
[(x-4)^2+y^2 <= 25, 9 <= x^2+(y-3)^2, 16 <= (x+w)^2+y^2],
[x, y], parameters = [w]);

which is called by procedure.

I hope that this bug will be addressed soon; then MultiIntPoly will work for the 32 bit version too.

@Avel Vogt:  are you on 32 bit too?

@max125 

plots:-animate(plot, [(x-1)*(x-2)*(x-3) + c, x=0.5 .. 3.5], c=0..1);

or

Explore(plot((x-1)*(x-2)*(x-3) + c, x=0.5 .. 3.5, view=-1..3), c=0. .. 1., animate);
(you may remove the animate option and use the slider instead).

@rlopez 

It is useful to note that if S₀ is a C^oo selection in an interval [a,b] then there is a continuous selection S in R which extends S₀; but a C^1 selection S extending S₀ may not exist. This is of course valid for eigenvalues too.

@digerdiga 

1. If we are interested only in real roots, the discontinuity cannot be avoided, as the example shows.

2. Yes, h(y)>0.384... is enough. The exact value is 2*sqrt(3)/9  but this is not important.

3. When all the complex roots are considered, the continuity can be obtained.
For example,
r1(c) = min({|z| : p(z) = 0},  r2(c) = max({|z| : p(z) = 0}
are continuous.

restart;
q:=(x-1)*(x-2)*(x-3):
sol:=c -> solve(q+c,x,explicit):
r1:=unapply( min(abs~([sol(c)])),c):
r2:=unapply( max(abs~([sol(c)])),c):
r:=c -> max(select(t->Im(t)=0, [sol(c)])):
plot(r, -8..8, discont, color=green);
plot([r1,r2], -8..8, color=[red,blue]);

@Rouben Rostamian  

I just want to mention that the piecewise method also handles non-convex domains (polygonal or not):

poly:= y<=x and y>=-x and (y>=4*x-3  or y<= -4*x+3) and x>=0 and x<=1;
F:=piecewise(poly, x^2+y^2, undefined);
G:=piecewise(poly, 3*sin(x)-y^2, undefined):
fieldplot([F,G], x =0 .. 1, y = -1 .. 1, arrows = SLIM, color = x) ;

@Carl Love 
Except that mine obtains the monomials of total degree exactly d.
To get "up to d" as OP wants it is even simpler:

mondeg:=proc(X::list(name),d)
expand((add(X)+1)^d);[op(%)]/~coeffs(%)
end:

Or, better (as you did):

mondegx:=proc(X::list(name),d)
local t;
coeffs(expand((add(X)+1)^d,X),X,t);[t]
end: