code: def num_jac(fun, t, y, f, threshold, factor, sparsity=None):
      """Finite differences Jacobian approximation tailored for ODE solvers.
  
      This function computes finite difference approximation to the Jacobian
      matrix of `fun` with respect to `y` using forward differences.
      The Jacobian matrix has shape (n, n) and its element (i, j) is equal to
      ``d f_i / d y_j``.
  
      A special feature of this function is the ability to correct the step
      size from iteration to iteration. The main idea is to keep the finite
      difference significantly separated from its round-off error which
      approximately equals ``EPS * np.abs(f)``. It reduces a possibility of a
      huge error and assures that the estimated derivative are reasonably close
      to the true values (i.e., the finite difference approximation is at least
      qualitatively reflects the structure of the true Jacobian).
  
      Parameters
      ----------
      fun : callable
          Right-hand side of the system implemented in a vectorized fashion.
      t : float
          Current time.
      y : ndarray, shape (n,)
          Current state.
      f : ndarray, shape (n,)
          Value of the right hand side at (t, y).
      threshold : float
          Threshold for `y` value used for computing the step size as
          ``factor * np.maximum(np.abs(y), threshold)``. Typically, the value of
          absolute tolerance (atol) for a solver should be passed as `threshold`.
      factor : ndarray with shape (n,) or None
          Factor to use for computing the step size. Pass None for the very
          evaluation, then use the value returned from this function.
      sparsity : tuple (structure, groups) or None
          Sparsity structure of the Jacobian, `structure` must be csc_matrix.
  
      Returns
      -------
      J : ndarray or csc_matrix, shape (n, n)
          Jacobian matrix.
      factor : ndarray, shape (n,)
          Suggested `factor` for the next evaluation.
      """
      y = np.asarray(y)
      n = y.shape[0]
      if n == 0:
          return np.empty((0, 0)), factor
  
      if factor is None:
          factor = np.full(n, EPS**0.5)
      else:
          factor = factor.copy()
  
      # Direct the step as ODE dictates, hoping that such a step won't lead to
      # a problematic region. For complex ODEs it makes sense to use the real
      # part of f as we use steps along real axis.
      f_sign = 2 * (np.real(f) >= 0).astype(float) - 1
      y_scale = f_sign * np.maximum(threshold, np.abs(y))
      h = (y + factor * y_scale) - y
  
      # Make sure that the step is not 0 to start with. Not likely it will be
      # executed often.
      for i in np.nonzero(h == 0)[0]:
          while h[i] == 0:
              factor[i] *= 10
              h[i] = (y[i] + factor[i] * y_scale[i]) - y[i]
  
      if sparsity is None:
          return _dense_num_jac(fun, t, y, f, h, factor, y_scale)
      else:
          structure, groups = sparsity
          return _sparse_num_jac(
              fun, t, y, f, h, factor, y_scale, structure, groups
          )
  
file_path: scipy/scipy/integrate/_ivp/common.py
imports: from itertools import groupby
  from warnings import warn
  import numpy
  from scipy.sparse import find, coo_matrix
label_desc: Use image node_id 7 for calling a global function with example usage: num_jac(fun, t, y, f, threshold, factor, sparsity) and returns: np,  factor, _dense_num_jac, _sparse_num_jac