flux integration ?

Dear fellows,

I'm trying to implement neumann boundary conditions on some edges on a domain of my CFD code we have this form of integration : http://www.uploadup.com/di-T10V.png where, "n" is normal vector and "Ni" is the shape function (Galerkin) now the question is how to solve this in MATLAB

lets say I'm trying to do the integration in X direction like: http://www.uploadup.com/di-FWPR.png