docs: def tri_diag_solve(a, b, c, d):
      """ Solve a[i] * x[i-1] + b[i] * x[i] + c[i] * x[i+1] = d[i] for x[i]
  
      Adapted from: http://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm.
      """
      n = len(a)
      a, b, c, d = map(numpy.array, (a, b, c, d))
      c[0] /= b[0]
      d[0] /= b[0]
      for i in range(1, n):
          m = b[i] - a[i] * c[i - 1]
          c[i] /= m
          d[i] = (d[i] - a[i] * d[i - 1]) / m
      x = d
      for i in range(n - 2, -1, -1):
          x[i] -= c[i] * x[i + 1]
      return x