code: def _complex_even(f, fx, x, h):
          
          n = len(x)
          ee = np.diag(h)
          hess = 2. * np.outer(h, h)
          for i in range(n):
              for j in range(i, n):
                  hess[i, j] = (f(x + 1j * ee[i] + ee[j])
                                - f(x + 1j * ee[i] - ee[j])).imag / hess[j, i]
                  hess[j, i] = hess[i, j]
          return hess
docs: Calculate Hessian with complex-step derivative approximation
  
          The stepsize is the same for the complex and the finite difference part