docs: def _solve_implicit_banded(current, banded_matrix):
      """Uses a banded solver for matrix inversion of a tridiagonal matrix.
  
      Converts the complete listed tridiagonal matrix *(nxn)* into a three row
      matrix *(3xn)* and calls :py:func:`scipy.linalg.solve_banded()`.
  
      :param array current:           the current state of the variable for which
                                      matrix inversion should be computed
      :param array banded_matrix:     complete diffusion matrix (*dimension: nxn*)
      :returns:                       output of :py:func:`scipy.linalg.solve_banded()`
      :rtype:                         array
  
      """
      #  can improve performance by storing the banded form once and not
      #  recalculating it...
      #  but whatever
      J = banded_matrix.shape[0]
      diag = np.zeros((3, J))
      diag[1, :] = np.diag(banded_matrix, k=0)
      diag[0, 1:] = np.diag(banded_matrix, k=1)
      diag[2, :-1] = np.diag(banded_matrix, k=-1)
      return solve_banded((1, 1), diag, current)