I have also noticed this issue from a yesterday post.

Download bugdiffint.mw

F := unapply(simplify(int(f, [yip = -infinity .. infinity, xip = -infinity .. infinity, tp = 0 .. t])),  t):
Optimization:-Maximize(F, 0..1/2);

I'll take a polynomial of degree 3 (it is easier to write here):
p(x) = x^3 + a1*x^2 + a2*x + a3.

The polynomial has 3 complex roots x1, x2, x3 which depend continuously on (a1,a2,a3) in C^3.
But this dependence is "global" i.e. the multi-function (a1,a2,a3) |--> {x1,x2,x3} is continuous.
If we want to obtain a continuous selection (as it is called in maths) we must jump from one branch to another.

Now, RootOf(p(x), x) denotes the first root (for evalf)  (see ?RootOf for the order of the roots)
and generally this is not continuous as a function of (a1,a2,a3).

A graphical illustration:

consider the polynomial p(x) = (x-1)*(x-2)*(x-3) + c  for 0 <= c <= 1
and look at its largest real root r(c); so, e.g. r(0)=3. The next animation shows clearly that r() is discontinuous as a function of c.

On the other side, if you have another (continuous) function h(y), there is the possibility that r(h(y)) be continuous [for example this happens here if h(y)>1].

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

For such problems it is much better to use a little math.

In this case, (1-x)^(10^9)  ≈ 1 - 10^9*x
with an excellent accuracy for a plot.

E1:=A+B+C+D:
E2:=A-B+C-D:
E3:=C:
sol:=solve(  [coeffs(expand(x*E1 + y*E2 + z*E3 - w*A), [A,B,C,D])], [x,y,z,w] )[]:
freevars:=indets(rhs~(sol),name):
[x,y,z,w]=~eval(rhs~(sol),freevars=~k);
# choose k for an integer solution (k=1 here); a simple procedure can be written for this if necessary

               [x = k, y = k, z = -2 k, w = 2 k]

Edit: Or, stealing Carl's isolve idea:

E1:=A+B+C+D:
E2:=A-B+C-D:
E3:=C:
sol:=isolve({coeffs(expand(x*E1 + y*E2 + z*E3 - w*A), [A,B,C,D])}):
freevars:=indets(rhs~(sol),name):
eval(sol, freevars =~1);
                 {w = 2, x = 1, y = 1, z = -2}

mondeg:=proc(X::list(name),d::posint)
if nops(X)=1 then X^~d else expand(add(X)^d); [op(%)]/~coeffs(%) fi
end:

This solution is longer but seems to be easier to follow.

Area:=proc(rel::set(relation))
local x,y, xy:=indets(rel,name),S, ff, u;
if nops(xy)<>2 then error "Two vars only!" fi;
x,y:=xy[];
S:=remove(hastype, [solve(rel)],`=`);
ff:=proc(rel)
local t,ux,uy,x1:=NULL,x2:=NULL,y1:=NULL,y2:=NULL, r:=rel;
ux:=select(t -> (lhs(t)=x and not has(rhs(t),y)) or (rhs(t)=x and not has(lhs(t),y)), rel);
uy:=select(t -> (lhs(t)=y and not has(rhs(t),x)) or (rhs(t)=y and not has(lhs(t),x)), rel);
if nops(uy)=2 and nops(ux)=0 then 
   ux,uy:=uy,ux; ux:=subs([x=y,y=x],ux); uy:=subs([x=y,y=x],uy); r:=subs([x=y,y=x],rel) fi;
if nops(ux)=2 and nops(uy)=0 then
   uy:=r minus ux;  if nops(uy)<>2 then error "Unusual solve solution" fi; 
   if rhs(ux[1])=x then x1:=lhs(ux[1]) fi;
   if rhs(ux[2])=x then x1:=lhs(ux[2]) fi;
   if lhs(ux[1])=x then x2:=rhs(ux[1]) fi;
   if lhs(ux[2])=x then x2:=rhs(ux[2]) fi;
   if rhs(uy[1])=y then y1:=lhs(uy[1]) fi;
   if rhs(uy[2])=y then y1:=lhs(uy[2]) fi;
   if lhs(uy[1])=y then y2:=rhs(uy[1]) fi;
   if lhs(uy[2])=y then y2:=rhs(uy[2]) fi;
   if nops([x1,x2,y1,y2])<4 then error "xk or yk null" fi;
   return simplify(int(1,y=y1..y2,x=x1..x2));
else error "Unable to compute!";
fi;   
end:
add(ff(u), u=S);
end:

R:={abs(abs(x-y)-abs(x+y)) >= 2*y-x+1, (x+1)^2+(y+1)^2 <= 2}:
Area(R);

radnormal(%);

evalf(%);
    4.419820813

[I could not upload the worksheet.]

It is a typesetting bug.

- If you execute
eval(%);

the strange term disappears.

- If you set
interface(typesetting=standard);
and re-execute, ==>

Error, (in print/diff) unable to display this derivative. Would you please report the problem to support@maplesoft.com

Yes, it's a bug.

I have renamed the variables (see the attach) and I wrote the system in a matrix form.

read "EQX.txt":
V:=[indets(EQX)[]]:nops(V);
                              1020
with(LinearAlgebra):
A,b:=LinearAlgebra:-GenerateMatrix( EQX, V ):
Rank(A);
                              326

Rank(<A|b>);
                              327

So, the system is incompatible (normally, Maple cannot be wrong for this computation, but we never know ...).

EQX.txt

The algorithm is simple.

Let {[x_k, u_k] :  k=1..n} union {[x_k, v_k] :  k=1..n}
be the extreme points given by convexhull.
Here u_k <= v_k, equality being possible.
(Obviously there cannot be 3 points with the same abscissa)
Also x_k are all distinct; one may suppose x_1 < ... < x_n
(after a sort command).

1. Compute min {u_k : k=1..n} = u_p
If p is not unique, take the largest one.

2. The wanted points are {[x_k, u_k] : k=p..n}.

You cannot use infolevel because the algorithms have changed. It seems that using another method (see ?Basis) does not help.
But the problem can be solved:

F:=[f[1],f[2],f[3]]:
gb:=Basis(F,plex(c,x[1],x[2])):  #  gbasis command is deprecated

factor~(gb);

From the last polynomial, it is clear that for c=3 we do not have a reduced basis.
So, if c is seen as a parameter, c=3 must be considered separately.

Y:=solve(sqrt(x^2+y^2)=(3-x)/2, y):
 X:=solve(Y[1]);
                           X := -3, 1
plots:-shadebetween(sqrt(x^2+y^2), (3-x)/2, x=X[1]..X[2], y=Y[1]..Y[2]);

Matrix(3, (i,j) -> Gradient~(V)[i][j]);

It will also work for non-cartesian coordinates.