Abstract Let M( f) denote the midpoint rule and T( f) the trapezoidal rule for estimating ∫ a b f(x) dx . Then Simpson's rule = λM( f)+(1− λ) T( f), where λ= 2 3 . We generalize Simpson's rule to multiple integrals as follows. Let D n be some polygonal region in R n, let P 0,…, P m denote the vertices of D n , and let P m+1 equal the center of mass of D n . Define the linear functionals M( f)=Vol( D n ) f( P n +1), which generalizes the midpoint rule, and T( f)=Vol( D n )([1/( m+1)]∑ j=0 m f( P j )), which generalizes the trapezoidal rule. Finally, our generalization of Simpson's rule is given by the cubature rule (CR) L λ = λM( f)+(1− λ) T( f), for fixed λ, 0⩽λ⩽1 . We choose λ, depending on D n , so that L λ is exact for polynomials of as large a degree as possible. In particular, we derive CRs for the n simplex and unit n cube. We also use points Q j ∈ ∂( D n ), other than the vertices P j , to generate T( f). This sometimes leads to better CRs for certain regions — in particular, for quadrilaterals in the plane. We use Grobner bases to solve the system of equations which yield the coordinates of the Q j 's.