Numerical double integrate in Java when boundaries are defined by x

Question

I am trying to calculate the integrate of the below expression

I tried to do it when Apache Commons Math library, but it does not solve double integrates.

 // c11, c12, c21 and c22 are constants.


 private double integrand(double x, double z) {
    return Math.exp(-0.5 *
            (matrix.getData()[0][0] * Math.pow(x, 2)) +
            (matrix.getData()[0][1] * z * x) +
            (matrix.getData()[1][0] * z * x) +
            (matrix.getData()[1][1] * Math.pow(z, 2))
    );
}

public double calculate() {
    int nx, ny;

    nx = (int) ((upperBoundX - lowerBoundX) / stepSizeX + 1);
    ny = (int) ((upperBoundY - lowerBoundY) / stepSizeY + 1);

    double[][] z = new double[nx][ny];
    double[] ax = new double[nx];
    double answer;

    for (int i = 0; i < nx; ++i) {
        for (int j = 0; j < ny; ++j) {
            z[i][j] = integrand(lowerBoundX + i * stepSizeX, lowerBoundY + j * stepSizeY);
        }
    }

    for (int i = 0; i < nx; ++i) {
        ax[i] = 0;
        for (int j = 0; j < ny; ++j) {
            if (j == 0 || j == ny - 1)
                ax[i] += z[i][j];
            else if (j % 2 == 0)
                ax[i] += 2 * z[i][j];
            else
                ax[i] += 4 * z[i][j];
        }
        ax[i] *= (stepSizeY / 3);
    }

    answer = 0;

    for (int i = 0; i < nx; ++i) {
        if (i == 0 || i == nx - 1) {
            answer += ax[i];
        }
        else if (i % 2 == 0) {
            answer += 2 * ax[i];
        }
        else {
            answer += 4 * ax[i];
        }
    }
    answer *= (stepSizeX / 3);

    return answer;
}

I tried to solve it with Simpson 1/3 rule but I could not make it work for integrates that have "x" in their boundaries. How can I solve this?