code: def linear_multistep_coeff(order, t, i, j, epsrel=1e-4):
      if order - 1 > i:
          raise ValueError(f"Order {order} too high for step {i}")
  
      def fn(tau):
          prod = 1.0
          for k in range(order):
              if j == k:
                  continue
              prod *= (tau - t[i - k]) / (t[i - j] - t[i - k])
          return prod
  
      return integrate.quad(fn, t[i], t[i + 1], epsrel=epsrel)[0]
  
imports: import torch
  from scipy import integrate
  from ...util import append_dims
label_desc: Use image node_id 1 for calling a global function with example usage: linear_multistep_coeff(order, t, i, j, epsrel) and returns: integrate, prod