code: def numerical_jacobian(func, input, eps=1e-3):
      output = func(input)
  
      ji = input.numpy().reshape(-1).shape[-1]
      jo = output.numpy().reshape(-1).shape[-1]
      NJ = np.zeros((jo, ji), dtype=np.float32)
  
      for i in range(ji):
          eps_perturb = mask_like(input.numpy(), i, mask_value=eps)
  
          output_perturb_add = (
              func(Tensor(input.numpy() + eps_perturb))
              .numpy()
              .reshape(-1)
          )
          output_perturb_sub = (
              func(Tensor(input.numpy() - eps_perturb))
              .numpy()
              .reshape(-1)
          )
  
          grad_approx = (
              (output_perturb_add) - (output_perturb_sub)
          ) / (2 * eps)
  
          NJ[:, i] = grad_approx
      return NJ
  
imports: import numpy
  from tinygrad.tensor import Tensor
label_desc: Use image node_id 3 for calling a global function with example usage: numerical_jacobian(func, input, eps) and returns: NJ