code: def _value_and_jacobian(func, xs, ys, xlims, ylims):
      """
      Compute *func* and its derivatives along x and y at positions *xs*, *ys*,
      while ensuring that finite difference calculations don't try to evaluate
      values outside of *xlims*, *ylims*.
      """
      eps = np.finfo(float).eps ** (
          1 / 2
      )  # see e.g. scipy.optimize.approx_fprime
      val = func(xs, ys)
      # Take the finite difference step in the direction where the bound is the
      # furthest; the step size is min of epsilon and distance to that bound.
      xlo, xhi = sorted(xlims)
      dxlo = xs - xlo
      dxhi = xhi - xs
      xeps = np.take([-1, 1], dxhi >= dxlo) * np.minimum(
          eps, np.maximum(dxlo, dxhi)
      )
      val_dx = func(xs + xeps, ys)
      ylo, yhi = sorted(ylims)
      dylo = ys - ylo
      dyhi = yhi - ys
      yeps = np.take([-1, 1], dyhi >= dylo) * np.minimum(
          eps, np.maximum(dylo, dyhi)
      )
      val_dy = func(xs, ys + yeps)
      return (val, (val_dx - val) / xeps, (val_dy - val) / yeps)
  
imports: import functools
  import numpy
  import matplotlib
  from matplotlib import _api
  from matplotlib.path import Path
  from matplotlib.transforms import Affine2D, IdentityTransform
  from .axislines import _FixedAxisArtistHelperBase, _FloatingAxisArtistHelperBase, GridHelperBase
  from .axis_artist import AxisArtist
  from .grid_finder import GridFinder
label_desc: Use image node_id 1 for calling a global function with example usage: _value_and_jacobian(func, xs, ys, xlims, ylims) and returns: val,  unknown,  unknown