def METHOD_NAME(fjac, X, D=None, **kwargs):
    nobs, ndim = X.shape
    if D is None:
        D = 0.1 * np.eye(ndim)
    elif np.isscalar(D):
        D = D * np.eye(ndim)
    U = np.zeros((nobs, 1))
    m = int(ndim * (ndim - 1) / 2)
    q0 = np.ones(m) * np.mean(np.diag(D)) * 1000
    for i in tqdm(range(nobs), "Calc Ao potential"):
        X_s = X[i, :]
        F = fjac(X_s)
        Q, _ = solveQ(D, F, q0=q0, **kwargs)
        H = np.linalg.inv(D + Q).dot(F)
        U[i] = -0.5 * X_s.dot(H).dot(X_s)
    return U.flatten(