Iterated Integrals

Code

graphicsFinal[sliceLocation_,showInnerIntegralPlot_,showSlice_]:=(
slice=If[showSlice,ParametricPlot3D[Append[xySlice[t,sliceLocation],s],{t,innerLeft,innerRight},{s,0,f@@xySlice[t,sliceLocation]},PlotStyle->Blue,Mesh->None,BoundaryStyle->Black],{}];
point=If[showInnerIntegralPlot,Graphics3D[{Blue,PointSize[0.02],Point[Append[xySlice[back,sliceLocation],NIntegrate[f@@xySlice[r,sliceLocation],{r,innerLeft,innerRight}]]]}],{}];
Show[innerIntegralPlot,innerIntegralCurve,point,surface,region,slice,Axes->Automatic,AxesLabel->{x,y,z},BoxRatios->{2,2,3},ViewPoint->{4,4,4},PlotRange->{{xBack,xRight},{yBack,yRight},{0,20}},Lighting->{{"Ambient",White}}]
)

graphicsInitial[xLeft_,xRight_,yLeft_,yRight_,innerIntegral_,showSurface_,showInnerIntegralPlot_,showSlice_,showRegion_]:=(

xBack=xLeft-(xRight-xLeft)/4;
yBack=yLeft-(yRight-yLeft)/4;
back=If[innerIntegral==1,xBack,yBack];
innerLeft=If[innerIntegral==1,xLeft,yLeft];
innerRight=If[innerIntegral==1,xRight,yRight];
outerLeft=If[innerIntegral==1,yLeft,xLeft];
outerRight=If[innerIntegral==1,yRight,xRight];

surface=If[showSurface,Plot3D[f[x,y],{x,xLeft,xRight},{y,yLeft,yRight},MeshStyle->Red,PlotStyle->Opacity[0],BoundaryStyle->Red],{}];
region=If[showRegion,Graphics3D[{Opacity[0],EdgeForm[Red],Polygon[{{xLeft,yLeft,0},{xLeft,yRight,0},{xRight,yRight,0},{xRight,yLeft,0}}]}],{}];
xySlice[inner_,outer_]:=Table[If[i==innerIntegral,inner,outer],{i,1,2,1}];

innerIntegralFunction[t_]:=Integrate[f@@xySlice[r,t],{r,innerLeft,innerRight}];

innerIntegralCurve=If[showInnerIntegralPlot,ParametricPlot3D[Append[xySlice[back,t],innerIntegralFunction[t]],{t,outerLeft,outerRight},PlotStyle->{Red,Thickness[0.008]}],{}];

innerIntegralPlot=If[showInnerIntegralPlot,ParametricPlot3D[Append[xySlice[back,t],s],{t,outerLeft,outerRight},{s,0,innerIntegralFunction[t]},PlotStyle->LightBlue,Mesh->None],{}];


)

graphics[xLeft_,xRight_,yLeft_,yRight_,innerIntegral_,sliceLocation_,showSurface_,showInnerIntegralPlot_,showSlice_,showRegion_]:=(
graphicsInitial[xLeft,xRight,yLeft,yRight,innerIntegral,showSurface,showInnerIntegralPlot,showSlice,showRegion];
graphicsFinal[sliceLocation,showSlice,showInnerIntegralPlot]
)

animation[xLeft_,xRight_,yLeft_,yRight_,innerIntegral_,showSurface_,showInnerIntegralPlot_,showSlice_,showRegion_,duration_,framerate_]:=(
graphicsInitial[xLeft,xRight,yLeft,yRight,innerIntegral,showSurface,showInnerIntegralPlot,showSlice,showRegion];
Table[Image[graphicsFinal[sliceLocation,showSlice,showInnerIntegralPlot]],{sliceLocation,outerLeft,outerRight,(outerRight-outerLeft)/(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.

f[x_,y_]:=3(x-1)^2+(y-2)^2+2;
xLeft=0;
xRight=2;
yLeft=0;
yRight= 3;

innerIntegral = 2;

sliceLocation=.5;
showSurface=True;
showInnerIntegralPlot=True;
showSlice = True;
showRegion=True;

graphics[xLeft,xRight,yLeft,yRight,innerIntegral,sliceLocation,showSurface,showInnerIntegralPlot,showSlice,showRegion]

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,yLeft,yRight,innerIntegral,showSurface,showInnerIntegralPlot,showSlice,showRegion,duration,framerate],DefaultDuration->duration]