code: def _multicomplex2(f, fx, x, h):
          
          n = len(x)
          ee = np.diag(h)
          hess = np.outer(h, h)
          cmplx_wrap = Bicomplex.__array_wrap__
          for i in range(n):
              for j in range(i, n):
                  zph = Bicomplex(x + 1j * ee[i, :], ee[j, :])
                  hess[i, j] = cmplx_wrap(f(zph)).imag12 / hess[j, i]
                  hess[j, i] = hess[i, j]
          return hess
docs: Calculate Hessian with Bicomplex-step derivative approximation