2D Simpson's Integrator

Fast, fully vectorised version of the Simpson's method for 2D domains. This code avoids the use of any for loops etc. For a given level of accuracy it can be an order of magnitude or more faster than dblquad.

Examples are provided within the code.

Useage:

ans = simp2D('func',xs,xe,ys,ye,NX,NY)

Input parameters:

func - 2D function that accepts vector inputs (otherwise
incorrect results can occur)
xs, xe - x-limits of integration
ys, ye - y-limits of integration
NX - number of integration intervals in the x-direction (should be even)
NY - number of integration intervals in the y-direction (should be even)