How to create a function that performs numerical double integration for a twovariable function?

I have tried to solve this assignment but I got stuck...

"Create a function that performs numerical double integration for a twovariable function. The input parameters are a handle of a function to be integrated and integration limits. The integration method is the rectangular method rectangular method. Compare performance of your function with the built-in function integral."

This is what I've got so far:

x1 = -10; x2 = 10;
f=@(x,y)x.^2+y.^2;
c=1;
n=500;
dx=(x2-x1)/n;
for i=0:dx:n-1
    y(c)=dx*f(1+i*dx)+dy*f(1+i*dy);
    c=c+1;
    s=sum(y)*dx;
    z=integral(f,x1,x2)
end

[Rik]

You're not making a function yet.

Also, since you're using the rectangular method, you can use ndgrid to generate a grid of coordinates within the range specified by the input. The interval size determines the weight of each value. So then you have all the data you need.

The call to the integral function should not be inside your function.

[per isakson]

See documentation on integral2, Numerically evaluate double integral