code: def convolve_functions(
      f, g, f_interval, g_interval, dt, fill=0, name=None, **kwargs
  ):
      r"""Expression containing numerical convolution of `fn1` with `fn2`
  
      Parameters
      ----------
      f : sympy expr
         An expression that is a function of t only.
      g : sympy expr
         An expression that is a function of t only.
      f_interval : (2,) sequence of float
         The start and end of the interval of t over which to convolve values of f
      g_interval : (2,) sequence of floats
         Start and end of the interval of t over which to convolve g
      dt : float
         Time step for discretization.  We use this for creating the
         interpolator to form the numerical implementation
      fill : None or float
         Value to return from sampling output `fg` function outside range.
      name : None or str, optional
         Name of the convolved function in the resulting expression.
         Defaults to one created by ``utils.interp``.
      \*\*kwargs : keyword args, optional
         Any other arguments to pass to the ``interp1d`` function in creating the
         numerical function for `fg`.
  
      Returns
      -------
      fg : sympy expr
         An symbolic expression that is a function of t only, and that can be
         lambdified to produce a function returning the convolved series from an
         input array.
  
      Examples
      --------
      >>> from nipy.algorithms.statistics.formula.formulae import Term
      >>> t = Term('t')
  
      This is a square wave on (0,1)
  
      >>> f1 = sympy.Piecewise((0, t <= 0), (1, t < 1), (0, True))
  
      The convolution of ``f1`` with itself is a triangular wave on [0, 2],
      peaking at 1 with height 1
  
      >>> tri = convolve_functions(f1, f1, [0, 2], [0, 2], 1.0e-3, name='conv')
  
      The result is a symbolic function
  
      >>> print(tri)
      conv(t)
  
      Get the numerical values for a time vector
  
      >>> ftri = lambdify(t, tri)
      >>> x = np.arange(0, 2, 0.2)
      >>> y = ftri(x)
  
      The peak is at 1
      >>> x[np.argmax(y)]
      1.0
      """
      # Note that - from the doctest above - y is
      """
      array([ -3.90255908e-16,   1.99000000e-01,   3.99000000e-01,
             5.99000000e-01,   7.99000000e-01,   9.99000000e-01,
             7.99000000e-01,   5.99000000e-01,   3.99000000e-01,
             1.99000000e-01,   6.74679706e-16])
      """
      # - so the peak value is 1-dt - rather than 1 - but we get the same
      # result from using np.convolve - see tests.
      f_vals = _eval_for(f, f_interval, dt)
      g_vals = _eval_for(g, g_interval, dt)
      fg_time, fg_vals = _conv_fx_gx(
          f_vals, g_vals, dt, min(f_interval), min(g_interval)
      )
      return interp(fg_time, fg_vals, fill=fill, name=name, **kwargs)
  
imports: import itertools
  import numpy
  import sympy
  from scipy.interpolate import interp1d
  from sympy import DiracDelta, Symbol
  from sympy.utilities.lambdify import implemented_function, lambdify
  from nipy.algorithms.statistics.formula.formulae import Formula, Term
label_desc: Use image node_id 11 for calling a global function with example usage: convolve_functions(f, g, f_interval, g_interval, dt, fill, name) and returns: interp