You can try if your version of MATLAB supports quad2d. This method uses adaptive simpsons quadrature to calculate 2D integrals.
Otherwise you need to refresh your knowledge in numerical methods. One solution may be to do as if you would calulate 2 integrals in one dimensions. You may define a grid and for this grid find all the function values for the grid. Then you integrate over the first dimension, let us call it x, for every grid point in y, similar to what you would have done by hand. For this you can use any intergral method available, eg. tarpetzoidal method or simpson quadrature. If the function is stiff, then you may want to use a runge kutta method.
Now you have a vector with the integral value of x for every point on y. This is analoge to having performed the integral analytically for f(x,y) over x to get a function g(y) = integral(f(x,y)) from a to b. Now, repeat the process in y. This time you only need to integrate once, and voila, you have a calculated the integral in 2 dimensions.
This link will show you a way to vectorize this operation:
Otherwise, a more intuitive code (which is slower than the vectorized would and thus should only be used to explain the theory and for programs where time is not an issue) would be
% Pseudo code
myGrid = myGridFun(inputs)
f = myFunctionValuesCalculation(myGrid);
g = zeros(size(f(1,:)));
for iter = 1:length(f(:,2))
g(iter) = myIntegrator(f(:,iter),binsize);
end
val = myIntegrator(g,binsize);
Where functions are used to make it easier to follow and is not a requirement.
As a final word I need to say that there probably are better methods than this for both speed and accuracy, but these would also require more work and knowledge. One of these methods is the Crank-Nicholson method that is a finite difference method used to solve partial differential equations. Given the grid size the method can of course be used on 2d integrals as well. The method described above is however a method that only require knowledge about 1 dimensional integrals and can be implemented quite easily on a few rows. Also, if access to any of MATLAB integration methods, these would most likely do better.
