diff --git a/doc/examples/schematics/L2.png b/doc/examples/schematics/L2.png
new file mode 100644
index 00000000..62597ee0
Binary files /dev/null and b/doc/examples/schematics/L2.png differ
diff --git a/doc/examples/schematics/L2.sch b/doc/examples/schematics/L2.sch
new file mode 100644
index 00000000..f738199d
--- /dev/null
+++ b/doc/examples/schematics/L2.sch
@@ -0,0 +1 @@
+L 1 2; right, inductors/coils=3, purple, thick
diff --git a/doc/tutorials.rst b/doc/tutorials.rst
index c042ec06..e0f95ee6 100644
--- a/doc/tutorials.rst
+++ b/doc/tutorials.rst
@@ -369,6 +369,7 @@ The input and output signals can be plotted using::
    >>> ax = v_i.plot((-1, 10), label='input')
    >>> ax = v_o.plot((-1, 10), axes=ax, label='output')
    >>> ax.legend()
+   >>> # Note, the show() method is not required when using IPython or Jupyter
    >>> plt.show()
 
 .. image:: examples/tutorials/basic/VRC2plot.png
diff --git a/lcapy/expr.py b/lcapy/expr.py
index b3d9f160..93ff3ed2 100644
--- a/lcapy/expr.py
+++ b/lcapy/expr.py
@@ -1701,6 +1701,9 @@ def convolve(self, x, commutate=False, **assumptions):
 
         If `commutate` is True, swap order of functions in integral.
 
+        `taumin` defaults to -oo and `taumax` defaults to oo.  These
+        are relaxed if either or both of the expressions are causal.
+
         The result is an unevaluated integral.  It can be evaluated using
         the `doit()` method.
 
diff --git a/lcapy/ltifilter.py b/lcapy/ltifilter.py
index 41834b61..f5905a1a 100644
--- a/lcapy/ltifilter.py
+++ b/lcapy/ltifilter.py
@@ -7,7 +7,7 @@
 from .expr import expr, equation, ExprTuple
 from .differentialequation import DifferentialEquation
 from .functions import Derivative, Function, exp, cos
-from .texpr import TimeDomainExpression
+from .texpr import texpr
 from .transfer import transfer
 from .symbols import t, s, omega0, j, pi, f
 from .utils import isiterable
@@ -177,7 +177,7 @@ def sdomain_initial_response(self, ic=None):
         from .sym import ssym
 
         if ic is None:
-            ic = ()
+            return s * 0
 
         if not isiterable(ic):
             ic = (ic, )
@@ -221,16 +221,19 @@ def initial_response(self, ic=None):
         Ysi = self.sdomain_initial_response(ic)
         return Ysi(t)
 
-    def response(self, x, ic=None, ni=None):
-        """Calculate response of filter to input `x` given a list of initial conditions
-        `ic` for time indexes specified by `ni`.  If `ni` is a tuple,
-        this specifies the first and last (inclusive) time index.
-
-        The initial conditions are valid prior to the time indices given by the ni
-        `x` can be an expression, a sequence, or a list/array of values.
+    def response(self, x, ic=None, use_time_domain=False):
+        """Calculate response of filter to input expression `x` given a list
+        of initial conditions `ic`.
         """
 
-        return 0
+        x = texpr(x)
+
+        if use_time_domain:
+            return self.impulse_response().convolve(x) + self.initial_response(ic)
+        else:
+
+            Y = self.transfer_function() * x(s) + self.sdomain_initial_response(ic)
+            return Y(t)
 
     def subs(self, *args, **kwargs):
 
diff --git a/lcapy/sexpr.py b/lcapy/sexpr.py
index 6f711e70..6c955922 100644
--- a/lcapy/sexpr.py
+++ b/lcapy/sexpr.py
@@ -493,7 +493,9 @@ def lti_filter(self, normalize_a0=True):
         return LTIFilter.from_transfer_function(self, normalize_a0)
 
     def dlti_filter(self, method='bilinear', alpha=0.5):
-        """Create DLTI filter using bilinear transform."""
+        """Create DLTI filter using bilinear transform.
+
+        See `discretize` for available methods."""
 
         from .transfer import transfer
 
@@ -764,6 +766,8 @@ def discretize(self, method=None, alpha=0.5, scale=None):
         'backward-diff', 'backward-euler' 'simpson', 'matched-Z',
         'zero-pole-matching'
 
+        See also `dlti_filter`.
+
         """
 
         signal = self.is_signal or self.is_squared or self.is_power