@robertocooper 
 

restart:
with(plots):
scheme1 := ["zgradient",["Blue","Cyan","Green","Yellow","Orange","Red"]]:
P1:=plot3d(x*y, colorscheme=scheme1,style=surfacecontour): 
P2:=plot3d(x*y,colorscheme=scheme1,style=point,symbol=asterisk): 
final_plot:=display({P1,P2},size=[2400,2400]);
plottools:-exportplot("abc2400.png",final_plot)

2400.pdf

@nm Unfortuntely the quality of the .eps and .pdf is very low (just zoom into your pdf). It is much better to use a high resolution .png.  It seems that the .eps is obtained by vectorizing a bitmap; this is the reason for its dimension (>3 MB) and quality.
The vector graphics obtained with Maxima are much better and smaller.

@minhthien2016 Hint. Using Pythagorean quadruples (see the wiki article) it's not very difficult to generate (all) such tetrahedrons for which the tangency points with the insphere have rational coordinates. It will suffice to take the center of the insphere in the origin. Try it.

@Carl Love I also don't know why you have presented the quantifiers here. The context seems to be clear: x, y are free variables. Maple is not designed to work with quantifiers. is and coulditbe are practically reserved for the assume facility.
Just try to formalize the definition of the uniform convergence in Maple (F(x,y) --> f(x), for y -->y0), using quantifiers.

@minhthien2016 Replace (in two places) ^+ with ^%T

(but then your version is older than 2018).

@Carl Love I think that it is simpler to consider 
diff( x = 1, x);

which produces the same 1=0.
But, is(x=1) and is(1=0)  are both false. That's OK, because  false implies false  simplifies to true. 

@minhthien2016 Please check the coonditions. Maybe there is a typo, but the method is obviously correct.

@Earl  Nice worksheet, congratulations.

@nm Yes, but the first version can e used e.g. in the header of a procedure. Not is also useful to construct other types such as And(symbol, Not(constant)).

@Earl In the animate command, the frames are constructed first, so a delay is not possible during the animation. We must generate more frames in the region where we want to slow down (just like in the procedure Q).

@Rouben Rostamian  I just want to add that this way we obtain just a particular solution, namely when u = u(x-mu*t).
It can be obtained with eval/subs:

eval(diff(u(x,t),t) +p*u(x,t)*diff(u(x,t),x) + q* diff(u(x,t),x$3), u(x,t)=U(x-mu*t)):
eval(%, x-mu*t=xi);  

For a general solution we should have used u(x,t) = U(xi,tau).

@Carl Love In this case, due to convexity, there will be either <= 6 solutions or infinitely many (if t=0).

@ruhamdam To obtain the polynomial combination, use

G, C := Basis(F, tdeg(t, u, v, w, x, y, z), G, output = extended)

I had not the patience to wait for the result ...
 

Edit. The order of the variables is of course essential. Compare e.g. quo(x^100+y, x+y, y)  with quo(x^100+y, x+y, x).

@ruhamdam OK, I see now. In my answer, it was not any difference between using the variables a[3],... and x3 ...
Now it is. The Groebner basis depends (sometimes drastically) on the monomial order, so it is important to know what is the substitution a[?]=x, ...; or maybe you changed the polynomials?

@ruhamdam In my answer, the second call to Basis was with the variables a[3],a[4],... substituted with x3, x4,... (without indices).

Your new worksheet also works for me in 1D; I did not try to find a formatting mistake in it, I have just converted it in 1D.
Please try the following:
Open a new worksheet with File > New > Worksheet Mode
and copy & paste the following:

with(Groebner);
with(PolynomialIdeals);
a := -9*w*v*t - 5/2*x*y*z*t + 6*w^3 + 5/2*x*z^2*w + 3*t^2*u + 3*y^3*t - 9*y*t - 4*y^2*z*w + 3/2*y*v*z^2 + 9*z*w - 1/2*z^3*u;
b := 2*x^4*z - 2*x^3*y^2 - 12*u^2*x*z + 36*u^2*y^2 - 48*u*v*x*y + 24*v^2*x^2 + 9*x^3 - 54*u^2;
c := 11*z^2*u*w + 20/3*x*z^2*y^2 - 12*t*y^2*v - 35*z*w^2*x + 24*y^2*z - 13/3*y^4*z - 90*w^2 - 11*u*y*z*t + 5*y*z*v*w + 10*x*v*t*z + 10*x*w*t*y - 7/3*x^3*x^2 + 22*w^2*y^2 + 3/2*t^2*x^2 - 15*x*z^2 - 3/2*v^2*z^2 + 36*t*v - 27*z;
d := 10*z*x*u*w + 10*y*u*v*z - 90*v^2 - 12*u*y^2*w + 20/3*z*x^2*y^2 + 11*x^2*v*t - 35*x*v^2*z + 22*v^2*y^2 + 3/2*u^2*z^2 - 7/3*x^3*z^2 - 3/2*x^2*w^2 + 36*u*w + 24*y^2*x - 13/3*y^4*x - 15*z*x^2 + 5*v*x*y*w - 11*x*y*u*t - 27*x;
e := 24*x*y*v*t - 87/2*x*v*z*w + 24*u*y*z*w + 33/2*x*u*z*t + 54*y^3 - 9*y^5 + 9/2*w*x^2*t + 51*y^2*v*w + 9/2*v*z^2*u - 7*x^2*z^2*y + 16*y^3*z*x - 45*u*y^2*t - 18*x*y*w^2 - 45*z*x*y - 18*y*v^2*z - 81*y + 81*u*t - 135*w*v;
f := 270*u^2*y*t - 3*x^2*y^2*w - 33/2*x^3*y*t + 57/2*x^3*z*w + 261*x*v^2*w - 72*u^2*z*w + 621*v*x*y + 54*x*u*w^2 + 405*z*x*u + 87/2*x^2*z^2*u - 219*y^3*v*x - 1134*u*y^2 + 270*u*y^4 + 108*x^2*w + 108*y*v^3 + 363/2*y*x^2*v*z - 171*x*u*v*t - 450*u*y*v*w - 285*u*y^2*z*x + 972*u;
g := 138*u*v*y*t - 129*x*y^2*v*z - 147*u*v*z*w + 187/2*x^2*y*z*w - 360*y^2*v + 62*v*y^4 + 27*x^2*t + 18*v^3*z - 48*x*u*w*t - 129/2*x*u*z^2*y - 243*u*y*z + 85*y^3*z*u - 21*x^2*y^2*t + 288*w*x*y + 288*z*x*v + 91/2*v*z^2*x^2 - 76*y^3*w*x - 24*x*v^2*t + 192*x*v*w^2 - 132*u*y*w^2 + 33*u^2*z*t + 9/2*x^3*z*t - 30*v^2*y*w + 486*v;
J := PolynomialIdeal(a, b, c, d, e, f, g);

then run it, pressing Enter. It should work.