In reply to Arne Henningsen.

Just to clarify, the quantity you are computing is a scalar, not a gradient. Regardless, you can use the QUAD subroutine in SAS/IML software to evaluate double integrals. See the article "Evaluate an iterated integral" and the SAS/IML documentation for the QUAD subroutine, which contains an example that is almost exactly what you need.