code: def newton_refine_solve(jac_both, x_val, surf_x, y_val, surf_y):
      r
      a_val, b_val, c_val, d_val = jac_both[:, 0]
      #       and
      e_val = x_val - surf_x
      f_val = y_val - surf_y
      # Now solve:
      denom = a_val * d_val - b_val * c_val
      delta_s = (d_val * e_val - c_val * f_val) / denom
      delta_t = (a_val * f_val - b_val * e_val) / denom
      return delta_s, delta_t
docs: r"""Helper for :func:`newton_refine`.
  
      We have a system:
  
      .. code-block:: rest
  
         [A C][ds] = [E]
         [B D][dt]   [F]
  
      This is not a typo, ``A->B->C->D`` matches the data in ``jac_both``.
      We solve directly rather than using a linear algebra utility:
  
      .. code-block:: rest
  
         ds = (D E - C F) / (A D - B C)
         dt = (A F - B E) / (A D - B C)
  
      Args:
          jac_both (numpy.ndarray): A ``4 x 1`` matrix of entries in a Jacobian.
          x_val (float): An ``x``-value we are trying to reach.
          surf_x (float): The actual ``x``-value we are currently at.
          y_val (float): An ``y``-value we are trying to reach.
          surf_y (float): The actual ``x``-value we are currently at.
  
      Returns:
          Tuple[float, float]: The pair of values the solve the
           linear system.