A brief animation to show the idea of iterated integration "in strips".

> map(with,{plottools,plots}):

>

> f:=(x,y)->(5/(x+1))*(2/(y+1));

> background:=plot3d(f(x,y),x=0..2,y=0..2,style=wireframe,
grid=[11,11], shading=ZHUE,axes=boxed):
background;

> dx:=.2:
dy:=.2:
xnum:=ceil(2/dx):
ynum:=ceil(2/dy):
xvals:=[seq(i*dx,i=0..xnum)]:
yvals:=[seq(i*dy,i=0..ynum)]:

> blocks:=[seq( seq( cuboid([xvals[i],yvals[j],0],
[xvals[i+1],yvals[j+1],f(0.5*dx+xvals[i],0.5*dy+yvals[j])]),
j=1..nops(yvals)-1),
i=1..nops(xvals)-1)]:

> frames:=[background,seq(display(blocks[1..3*k],background,
view=[0..2,0..2,0..10]),k=1..nops(blocks)/3),
display(blocks,background,view=[0..2,0..2,0..10])]:

> display(frames,insequence=true,view=[0..2,0..2,0..10],
title="y first int",labels=["x","y","z"]);

> cubevol:=add(add(f(0.5*dx+xvals[i],0.5*dy+yvals[j])*dx*dy,
j=1..nops(yvals)-1),i=1..nops(xvals)-1);

> realvol:=evalf(int(int(f(x,y),y=0..2),x=0..2));

>