Integrating Over General Regions

Code

graphicsFinal[sliceLocation_,showSlice_]:=(
bottomSolutions=NSolve[{Rationalize[bottomBound[x]==sliceLocation],xLeft≤x≤xRight},x];
topSolutions=NSolve[{Rationalize[topBound[x]==sliceLocation],xLeft≤x≤xRight},x];
leftBound=If[innerIntegral==1,If[topSolutions=={},left,x/.First[topSolutions]],bottomBound[sliceLocation]];
rightBound=If[innerIntegral==1,If[bottomSolutions=={},right,x/.First[bottomSolutions]],topBound[sliceLocation]];
slice=If[showSlice,ParametricPlot3D[Append[xySlice[t,sliceLocation],s],{t,leftBound,rightBound},{s,0,f@@xySlice[t,sliceLocation]},PlotStyle->Blue,Mesh->None,BoundaryStyle->Black],{}];

Show[surface,region,slice,surfaceAboveRegion,Axes->Automatic,AxesLabel->{x,y,z},BoxRatios->{2,2,3},ViewPoint->{-4,4,4},PlotRange->{{xLeft,xRight},{yLeft,yRight},{0,zMax}},Lighting->{{"Ambient",White}}]
)

graphicsInitial[xLeft_,xRight_,innerIntegral_,showEntireSurface_,showRegionOfIntegration_,showSurfaceAboveRegionOfIntegration_,showSlice_]:=(

yLeft=NMinimize[{bottomBound[x],xLeft≤x≤xRight},x][[1]];
yRight=NMaximize[{topBound[x],xLeft≤x≤xRight},x][[1]];

left=If[innerIntegral==1,yLeft,xLeft];right=If[innerIntegral==1,yRight,xRight];

zMax=NMaximize[{f[x,y],xLeft≤x≤xRight,yLeft≤y≤yRight},{x,y}][[1]];

surface=If[showEntireSurface,Plot3D[f[x,y],{x,xLeft,xRight},{y,yLeft,yRight},MeshStyle->{{Red,Opacity[0.2]}},PlotStyle->Opacity[0],BoundaryStyle->None],{}];
region=If[showRegionOfIntegration,ParametricPlot3D[{s,t,0},{s,xLeft,xRight},{t,bottomBound[s],topBound[s]},PlotStyle->LightGray,Mesh->None,BoundaryStyle->{Black,Thickness[0.005]}],{}];
surfaceAboveRegion=If[showSurfaceAboveRegionOfIntegration,ParametricPlot3D[{s,t,f[s,t]},{s,xLeft,xRight},{t,bottomBound[s],topBound[s]},MeshStyle->{{Red,Thickness[0.005]}},PlotStyle->Opacity[0],BoundaryStyle->{Red,Thickness[0.005]}],{}];
xySlice[inner_,outer_]:=Table[If[i==innerIntegral,inner,outer],{i,1,2,1}];
)

graphics[xLeft_,xRight_,innerIntegral_ ,sliceLocation_,showEntireSurface_,showRegionOfIntegration_,showSurfaceAboveRegionOfIntegration_,showSlice_]:=(
graphicsInitial[xLeft,xRight,innerIntegral,showEntireSurface,showRegionOfIntegration,showSurfaceAboveRegionOfIntegration,showSlice];
graphicsFinal[sliceLocation,showSlice]
)

animation[xLeft_,xRight_,innerIntegral_,showEntireSurface_,showRegionOfIntegration_,showSurfaceAboveRegionOfIntegration_,showSlice_,duration_,framerate_]:=(
graphicsInitial[xLeft,xRight,innerIntegral,showEntireSurface,showRegionOfIntegration,showSurfaceAboveRegionOfIntegration,showSlice];
epsilon=(right-left)/1000;
Table[Image[graphicsFinal[sliceLocation,showSlice]],{sliceLocation,left+epsilon,right-epsilon,(right-left)/(duration*framerate)}]
)

Graphics

Modify the functions and options below to change the graphics.
innerIntegral = 1 specifies that the inner integral of the iterated integrals should be with respect to x, and innerIntegral = 2 specifies that the inner integral is with respect to y.
Warning: The code makes some assumptions about the behavior of bottomBound and topBound.  In order to avoid errors, bottomBound and topBound must be functions that do not intersect between the endpoints.  Additionally, if using innerIntegral=1, they should both be increasing on the interval [xLeft, xRight].

f[x_,y_]:=x^2 y;
bottomBound[x_]:=x^3;
topBound[x_]:=Sqrt[x];

xLeft=0;
xRight=1;

innerIntegral = 2;

sliceLocation=0.5;

showEntireSurface=True;
showRegionOfIntegration=True;
showSurfaceAboveRegionOfIntegration=True;
showSlice=True;

graphics[xLeft,xRight,innerIntegral,sliceLocation,showEntireSurface,showRegionOfIntegration,showSurfaceAboveRegionOfIntegration,showSlice]

Animation

All options except sliceLocation are taken from the Graphics section above.
It may take some time to run this command, as all frames of the animation are precomputed so that it can run smoothly.

duration=3;
framerate=16;

ListAnimate[animation[xLeft,xRight,innerIntegral,showEntireSurface,showRegionOfIntegration,showSurfaceAboveRegionOfIntegration,showSlice,duration,framerate],DefaultDuration->duration]