By coloring the cube with color f(x,y,z) you will have information on f only on the boundary of the cube.
That is, only f(0,y,z), f(1,y,z),...,f(x,y,1) are visualized.
You can't "see" e,g. the value f(0.5, 0.6, 0.7).
So, you will see the max/min only if they are on the boundary.

@Carl Love 

Sorry for the confusion, but I thought that the reference to Eigenvectors is obvious because of the presence of lambda.

@Carl Love 

Eigenvectors produces a random order for the eigenvalues because after solving the characteristic equation with the local variable lambda, indices is called.
You are right about indexorder. The procedure which calls indices without indexorder is `convert/multiset`. If we change it inserting indexorder, the random behavior disappears.

@Markiyan Hirnyk 

OK, it's enough for me too.
P.S. I wish you good luck in obtaining small Groebner bases using the cited articles.

@Markiyan Hirnyk 

I was not saying that the paper is not serious. But their approach is only for (some) numeric problems. The Groebner bases are used for much more!

@Markiyan Hirnyk 

I did not read them because as I have added (probably you did not see because it was during your reply) I am interested in general Groebner bases which are inherently huge.

Edit. OK, you made me curious and I have looked over [1]. As I have anticipated, they approximate a Groebner basis using Newton's method. That is not what most people want!

@Markiyan Hirnyk 

OK, but if for 70 equations we need supercomputers then ...

Edit. I was (and am) interested in general Groebner bases, not in some versions designed for numerics. Here the situation is unfortunately clear: the Grobner basis is unique and very often is HUGE.

@Markiyan Hirnyk 

The Groebner basis technique is excellent in theory and works for small size problems. I like it very much.
Unfortunately for even moderately large problems it cannot be used. The compexity could be double exponential O(a^(b^n))   !!

@Carl Love 

Yes, irem is a little faster than mods, but if p and nops(S) are huge,
taking (at least from time to time) the irem of the partial product
will speed up the computations.

@acer 
Probably a more instructive example would be:

progbar:= proc(p::realcons)
local x;
if p>100 or p<0 then DocumentTools:-Tabulate( [[]] ): return NULL fi:
DocumentTools:-Tabulate([ sprintf("%a %c", p, "%"),
                          plot(100,x=0..p, view=[0..100,0..100],
                                  filled,axes=box, tickmarks=[[],[]],labels=["",""], size=[450,50]) ],
                          exterior=none, interior=none, widthmode=pixels, width=450 );
end:

s:=0;
"starting ...";
for i to 10000 do
  s:=s + evalf(1/i);
  if i mod 100 = 0 then  progbar(i/100) fi;
od:
progbar(101):
's'=s;

@dpaddy 

If g is continuous and C^1 piecewise continuous then you can use IntegrationTools[Split].
If g is given as piecewise, Maple should be able do compute directly int(f*D(g), a..b);

If g is C^1 piecewise continuous but not (globally) continuous you must also consider the jumps.

@dpaddy 

If in the integral of f*dg  one of the functions  f,g  is C^1 (at least piecewice) then the integral can be reduced to a Riemann one
and it is not difficult to write a Maple program to do this.
(But be aware of the fact that there exists (pathological) functions for which Maple fails in computing the Riemann integral.)

Otherwise, you have to study theoretically the problem, and if the integral exists try to reduce it to the previous case (if possible).
Note that the Stieltjes integral is not very simple in the general case. Even if both f,g are continuous, it may not exist.

@AmusingYeti 

Any 3d plot!

plot3d(x^2-y^2,x=-1..1, y=-1..1);
#Here is the print screen; I have posted such images here long time ago.
(I have to  use Maxima for vector 3d plots!)

@AmusingYeti 

It must be added that unfortunately the .eps export for 3d graphics is terrible (unusable), at least in Windows.

@Carl Love 

g(x)=floor(x):

inf(x*dg(x), x=0..5) =

x*g(x)|(x=5) - x*g(x)|(x=0)  - int(g(x), x=0..5) =

25 - 10 =

15

Now, using Maple:

Int(x*D(g)(x), x=0..5):
IntegrationTools[Parts](%,x):
eval(%, g=floor):
value(%);
                               15

# Wrong result if trying directly!
int(x*D(floor)(x), x=0..5);
                               0