code: def _integrate_rk12(
      x0, y0, dmap, f, maxlength, broken_streamlines=True
  ):
      """
      2nd-order Runge-Kutta algorithm with adaptive step size.
  
      This method is also referred to as the improved Euler's method, or Heun's
      method. This method is favored over higher-order methods because:
  
      1. To get decent looking trajectories and to sample every mask cell
         on the trajectory we need a small timestep, so a lower order
         solver doesn't hurt us unless the data is *very* high resolution.
         In fact, for cases where the user inputs
         data smaller or of similar grid size to the mask grid, the higher
         order corrections are negligible because of the very fast linear
         interpolation used in `interpgrid`.
  
      2. For high resolution input data (i.e. beyond the mask
         resolution), we must reduce the timestep. Therefore, an adaptive
         timestep is more suited to the problem as this would be very hard
         to judge automatically otherwise.
  
      This integrator is about 1.5 - 2x as fast as RK4 and RK45 solvers (using
      similar Python implementations) in most setups.
      """
      # This error is below that needed to match the RK4 integrator. It
      # is set for visual reasons -- too low and corners start
      # appearing ugly and jagged. Can be tuned.
      maxerror = 0.003
  
      # This limit is important (for all integrators) to avoid the
      # trajectory skipping some mask cells. We could relax this
      # condition if we use the code which is commented out below to
      # increment the location gradually. However, due to the efficient
      # nature of the interpolation, this doesn't boost speed by much
      # for quite a bit of complexity.
      maxds = min(1.0 / dmap.mask.nx, 1.0 / dmap.mask.ny, 0.1)
  
      ds = maxds
      stotal = 0
      xi = x0
      yi = y0
      xyf_traj = []
  
      while True:
          try:
              if dmap.grid.within_grid(xi, yi):
                  xyf_traj.append((xi, yi))
              else:
                  raise OutOfBounds
  
              # Compute the two intermediate gradients.
              # f should raise OutOfBounds if the locations given are
              # outside the grid.
              k1x, k1y = f(xi, yi)
              k2x, k2y = f(xi + ds * k1x, yi + ds * k1y)
  
          except OutOfBounds:
              # Out of the domain during this step.
              # Take an Euler step to the boundary to improve neatness
              # unless the trajectory is currently empty.
              if xyf_traj:
                  ds, xyf_traj = _euler_step(xyf_traj, dmap, f)
                  stotal += ds
              break
          except TerminateTrajectory:
              break
  
          dx1 = ds * k1x
          dy1 = ds * k1y
          dx2 = ds * 0.5 * (k1x + k2x)
          dy2 = ds * 0.5 * (k1y + k2y)
  
          ny, nx = dmap.grid.shape
          # Error is normalized to the axes coordinates
          error = np.hypot(
              (dx2 - dx1) / (nx - 1), (dy2 - dy1) / (ny - 1)
          )
  
          # Only save step if within error tolerance
          if error < maxerror:
              xi += dx2
              yi += dy2
              try:
                  dmap.update_trajectory(xi, yi, broken_streamlines)
              except InvalidIndexError:
                  break
              if stotal + ds > maxlength:
                  break
              stotal += ds
  
          # recalculate stepsize based on step error
          if error == 0:
              ds = maxds
          else:
              ds = min(maxds, 0.85 * ds * (maxerror / error) ** 0.5)
  
      return stotal, xyf_traj
  
imports: import numpy
  import matplotlib
  from matplotlib import _api, cm, patches
  import matplotlib.colors
  import matplotlib.collections
  import matplotlib.lines
label_desc: Use image node_id 3 for calling a global function with example usage: _integrate_rk12(x0, y0, dmap, f, maxlength, broken_streamlines) and returns: stotal,  xyf_traj