2D Simpsons Composite integration

In the Mathematics Stack Exchange I questioned about the extension of the composite Simpsons rule to 2D integration.

The problem is that Simpsons rule requires and even number of subintervals and I'm concerned with the possibility of having an even number of subintervals. In another forum, for 1D integration with Simpsons rule, the suggestion was to take out the last points and use Simpsons 3/8 rule on this last 4 points to keep the integration in the same level of accuracy (the same error order of .

My question on the forum was about the extension of this approach to a 2D integration, if the weights were corrected. Assuming they are correct, I coded the function:

function I = simpson2D(x,y,z)
% Composite Simpson's rule for 2D integration

    from meshgrid
    Nx = size(x,2); % number of columns
    Ny = size(x,1); % number of rows

    hx = x(1,2)-x(1,1);
    hy = y(2,1)-y(1,1);

    if (mod(Nx,2)) % Nx is odd (even number of intervals)
        NxIsEven = false;
        nx = Nx;
        wx = ones(1,Nx);
        wx(2:2:nx-1) = 4;
        wx(3:2:nx-2) = 2;
    else % Nx is even (odd number of intervals)
        NxIsEven = true;
        nx = Nx-3;
        wx = ones(1,Nx);
        wx(2:2:nx-1) = 4;
        wx(3:2:nx-2) = 2;
        wx(Nx-2:Nx) = [3 3 1];
    end

    if (mod(Ny,2)) % Ny is odd (even number of intervals)
        NyIsEven = false;
        ny = Ny;
        wy = ones(Ny,1);
        wy(2:2:ny-1) = 4;
        wy(3:2:ny-2) = 2;
    else % Ny is even (odd number of intervals)
        NyIsEven = true;
        ny = Ny-3;
        wy = ones(Ny,1);
        wy(2:2:ny-1) = 4;
        wy(3:2:ny-2) = 2;
        wy(Ny-2:Ny) = [3 3 1];
    end

    w = wy*wx;

    if (NxIsEven && NyIsEven)
        I = sum(sum(w(1:ny,1:nx).*z(1:ny,1:nx)))/9 + ...
           (sum(sum(w(1:ny,nx:Nx).*z(1:ny,nx:Nx)))+sum(sum(w(ny:Ny,1:nx).*z(ny:Ny,1:nx))))/8 + ...
            sum(sum(w(ny:Ny,nx:Nx).*z(ny:Ny,nx:Nx)))*9/64;
    elseif (~NxIsEven && NyIsEven)
        I = sum(sum(w(1:ny,1:nx).*z(1:ny,1:nx)))/9 + ...
            sum(sum(w(ny:Ny,1:nx).*z(ny:Ny,1:nx)))/8;
    elseif (NxIsEven && ~NyIsEven)
        I = sum(sum(w(1:ny,1:nx).*z(1:ny,1:nx)))/9 + ...
            sum(sum(w(1:ny,nx:Nx).*z(1:ny,nx:Nx)))/8;
    else % (~NxIsEven && ~NyIsEven)
        I = sum(sum(w(1:ny,1:nx).*z(1:ny,1:nx)))/9;
    end
    I = I*hx*hy;
end

The question is simple: can I improve this code?

How could I write it without so many if-elseif to account for the different possibilities on the parity of Nx and Ny? How to better define the weights vectors wx and inwy?