public static void hessianCD(final ObjectiveFunctionFast afunc, final double[] ax, final double afx, double[][] ahess) {
        int n = afunc.getDim();
        double[] h = new double[n];
        double xh, oldx, oldy;
        int i, j;
        double tol = Math.pow(MachineAccuracy.EPSILON, 0.25);//TODO revise maybe 1/3
        double[] fplus = new double[n];
        double[] fminus = new double[n];

        double fxx;

        for (i = 0; i < n; i++) {
            h[i] = tol * (Math.abs(ax[i]) + MachineAccuracy.EPSILON);
            xh = ax[i] + h[i];
            h[i] = xh - ax[i];
            oldx = ax[i];
            ax[i] = oldx + h[i];
            fplus[i] = afunc.f(ax);
            ax[i] = oldx - h[i];
            fminus[i] = afunc.f(ax);
            ax[i] = oldx;
        }

        double tH;
        for (i = 0; i < n; i++) {
            ahess[i][i] = h[i] * h[i];
            for (j = i + 1; j < n; j++) {
                tH = h[i] * h[j];
                ahess[i][j] = tH;
                ahess[j][i] = tH;
            }
        }

        for (i = 0; i < n; i++) {
            for (j = 0; j < n; j++) {
                if (i == j) {
                    ahess[i][j] = (fplus[i] + fminus[j] - 2 * afx) / ahess[i][j];
                } else {
                    oldx = ax[i];
                    oldy = ax[j];
                    ax[i] = oldx + h[i];
                    ax[j] = oldy - h[j];
                    fxx = afunc.f(ax);
                    ax[i] = oldx;
                    ax[j] = oldy;
                    ahess[i][j] = (fplus[i] + fminus[j] - afx - fxx) / ahess[i][j];
                }
            }
        }
        //H=(H+H')/2
        for (i = 0; i < n; i++) {
            for (j = i + 1; j < n; j++) {
                tH = (ahess[i][j] + ahess[j][i]) / 2.0;
                ahess[i][j] = tH;
                ahess[j][i] = tH;
            }
        }
    }