> #The region files are graphs of regions used in triple integration
> #This one is for the region -1<=x<=1,-1<=y<=1,-1<=z<=1
> #The general form for x: a<=x<=b (a is left, b is right)
> # for y: f(x)<=y<=g(x) (f(x) is front, g(x) is rear) 
> # for z: F(x,y)<=z<=G(x,y) (F(x,y) is bottom, G(x,y) is top)
> a:=-1;b:=1;#x limits x=a to x=b
> f:=-1;g:=1;#y limits y=f(x) to y=g(x)
> F:=-1;G:=1;#z limits z =F(x,y) to z=G(x,y)
> triple:=Int(Int(Int(1,z=F..G),y=f..g),x=a..b);
> n:=21;m:=9;#n is the number of x divisions, m the number in the y divsions
> xstep:=a+i*(b-a)/n:x0:=eval(xstep,i=0):xn:=eval(xstep,i=n):
> ystep:=eval(f+j*(g-f)/m,{x=xstep}):y0:=eval(ystep,j=0):ym:=eval(ystep,j=m):ystep0:=eval(ystep,i=0):ystepn:=eval(ystep,i=n):
> Top:=[seq([seq(eval([xstep, ystep, G],{x=xstep,y=ystep}),j=0..m)],i=0..n)]:
> Bottom:=[seq([seq(eval([xstep, ystep,F],{x=xstep,y=ystep}),j=0..m)],i=0..n)]:
> Front:=[seq(eval([[xstep,y0,G],[xstep,y0,F]],{x=xstep,y=y0}),i=0..n)]:
> Rear:=[seq(eval([[xstep,ym,G],[xstep,ym,F]],{x=xstep,y=ym}),i=0..n)]:
> Left:=[[seq(eval([x0,ystep,F],{x=x0,y=ystep0,i=0}),j=0..m)], [seq(eval([x0, ystep, G],{x=x0,y=ystep0,i=0}),j=0..m)]]:
> Right:=[[seq(eval([xn,ystep,F],{x=xn,y=ystepn,i=n}),j=0..m)], [seq(eval([xstep, ystep, G],{x=xn,y=ystepn,i=n}),j=0..m)]]:
> PLOT3D(MESH(evalf(Top)),MESH(evalf(Bottom)),MESH(evalf(Front)),MESH(evalf(Rear)),MESH(Left),MESH(Right));
>