code: def test_num_jac():
      def fun(t, y):
          return np.vstack(
              [
                  -0.04 * y[0] + 1e4 * y[1] * y[2],
                  0.04 * y[0] - 1e4 * y[1] * y[2] - 3e7 * y[1] ** 2,
                  3e7 * y[1] ** 2,
              ]
          )
  
      def jac(t, y):
          return np.array(
              [
                  [-0.04, 1e4 * y[2], 1e4 * y[1]],
                  [0.04, -1e4 * y[2] - 6e7 * y[1], -1e4 * y[1]],
                  [0, 6e7 * y[1], 0],
              ]
          )
  
      t = 1
      y = np.array([1, 0, 0])
      J_true = jac(t, y)
      threshold = 1e-5
      f = fun(t, y).ravel()
  
      J_num, factor = num_jac(fun, t, y, f, threshold, None)
      assert_allclose(J_num, J_true, rtol=1e-5, atol=1e-5)
  
      J_num, factor = num_jac(fun, t, y, f, threshold, factor)
      assert_allclose(J_num, J_true, rtol=1e-5, atol=1e-5)
  
imports: from itertools import product
  from numpy.testing import assert_, assert_allclose, assert_array_less, assert_equal, assert_no_warnings, suppress_warnings
  import pytest
  from pytest import raises
  import numpy
  from scipy.optimize._numdiff import group_columns
  from scipy.integrate import solve_ivp, RK23, RK45, DOP853, Radau, BDF, LSODA
  from scipy.integrate import OdeSolution
  from scipy.integrate._ivp.common import num_jac
  from scipy.integrate._ivp.base import ConstantDenseOutput
  from scipy.sparse import coo_matrix, csc_matrix
label_desc: Use image node_id 38 for calling a global function with example usage: test_num_jac() and returns: np, np