A = Int(f(x),x=a..b); f := (x,y) -> 3*x^2+(y-2)^2+2: vf := f(x,y): 'f(x,y)' = f(x,y); Int('f',A) = Int( Int('f(x,y)',x=0..1),y=0..3); Int( Int('f(x,y)',x=0..1),y=0..3) = int( int(f(x,y),x=0..1),y=0..3); plot3d(vf,x=0..1,y=0..3,axes=frame); Rsum1 := f(1/2,1/2)+f(1/2,3/2)+f(1/2,5/2); ApproxInt2d := proc(f,xexpress,yexpress,n,m) local a,b,c,d,dx,i,j,dy; a := op(1,op(2,xexpress)); b := op(2,op(2,xexpress)); c := op(1,op(2,yexpress)); d := op(2,op(2,yexpress)); dx := evalf((b-a)/n); dy := evalf((d-c)/m); sum(sum(subs(x = a+(i-.5)*dx,y = c+(j-.5)*dy,f),i = 1 .. n),j = 1 .. m)*dx*dy end: Rsum1_Now := ApproxInt2d(f(x,y),x=0..1,y=0..3,1,3); Rsum2 := ApproxInt2d(f(x,y),x=0..1,y=0..3,2,3); RS := proc(m,n) ApproxInt2d(f(x,y),x=0..1,y=0..3,m,n) end: RS(1,3), RS(2,3), RS(3,4), RS(10,10), RS(100,100); RS(1,3), RS(1,30), RS(1,60), RS(1,200), RS(1,2000); Int('f(1/2,y)',y=0..3) = Int(f(1/2,y),y=0..3); Int(f(1/2,y),y=0..3) = int(f(1/2,y),y=0..3); evalf(45/4); Int(f(x,y),y=0..3) = int(f(x,y),y=0..3); Int(9*x^2+9,x=0..1) = int(9*x^2+9,x=0..1); #