How can I evaluate the double integration of the stiffness matrix using Gaussian quadrature in Matlab?

Hi there everybody. In finite elements analysis, the stiffness matrix using Galerkin method is found to be expressed as this: K(i,j)=int(int(-D*((d2Ndx2(j)*d2Ndx2(i))+ v*(d2Ndy2(j)*d2Ndx2(i))+(d2Ndy2(j)*d2Ndy2(i))+ v*(d2Ndx2(j)*d2Ndy2(i))+(2*(1-v)*(d2Ndxdy(j)*d2Ndxdy(i))))dx)dy). N is a vector containing the shape function D and v are material properties. I'm using 4 node rectangular element with 3 degrees of freedom at each node which means ( N ) is a vector of 12 shape functions. Can anyone help me writing a Matlab code to evaluate the stiffness matrix using Gaussian quadrature rule or at least show me how to transform the second derivatives from (x-y)domain to (xi-eta)domain . thanks a lot.