diff --git a/kwant/continuum/_common.py b/kwant/continuum/_common.py
index efa5b1260ed8fa02fcb5dd9b7c2adfe9dbe00530..58bc4b6dae79ffddff1d35c8af8046b1262a01e8 100644
--- a/kwant/continuum/_common.py
+++ b/kwant/continuum/_common.py
@@ -39,38 +39,35 @@ del default_subs['I']
 del default_subs['pi']
 
 
-################  Various helpers to handle sympy
+################  Helpers to handle sympy
 
+def lambdify(expr, subs=None):
+    """Return a callable object for computing a continuum Hamiltonian.
 
-def lambdify(hamiltonian, subs=None):
-    """Return a callable object for computing continuum Hamiltonian.
+    .. warning::
+        This function uses ``eval`` (because it calls ``sympy.sympify``), and
+        thus should not be used on unsanitized input.
 
     Parameters
     ----------
-    hamiltonian : sympy.Expr or sympy.Matrix, or string
-        Symbolic representation of a continous Hamiltonian. When providing
-        a sympy expression, it is recommended to use ``momentum_operators``
-        defined in `~kwant.continuum` in order to provide
-        proper commutation properties. If a string is provided it will be
-        converted to a sympy expression using `kwant.continuum.sympify`.
-    subs: dict, defaults to empty
-        A namespace of substitutions to be performed on the input
-        ``hamiltonian``. It can be used to simplify input of matrices or
-        alternate input before proceeding further. Please see examples below.
+    expr : str or SymPy expression
+        Expression to be converted into a callable object
+    subs : dict or ``None`` (default)
+        Namespace of substitutions to be performed on `expr`.
 
     Example:
     --------
-        >>> f = kwant.continuum.lambdify('a + b', subs={'b': 'b + c'})
+        >>> f = lambdify('a + b', subs={'b': 'b + c'})
         >>> f(1, 3, 5)
         9
 
-        >>> subs = {'s_z': [[1, 0], [0, -1]]}
-        >>> f = kwant.continuum.lambdify('k_z**2 * s_z', subs=subs)
-        >>> f(.25)
-        array([[ 0.0625,  0.    ],
-               [ 0.    , -0.0625]])
+        >>> subs = {'sigma_plus': [[0, 2], [0, 0]]}
+        >>> f = lambdify('k_x**2 * sigma_plus', subs)
+        >>> f(0.25)
+        array([[ 0.   ,  0.125],
+               [ 0.   ,  0.   ]])
     """
-    expr = sympify(hamiltonian, subs)
+    expr = sympify(expr, subs)
 
     args = [s.name for s in expr.atoms(sympy.Symbol)]
     args += [str(f.func) for f in expr.atoms(AppliedUndef, sympy.Function)]
@@ -97,20 +94,27 @@ def sympify(expr, subs=None):
     In addition, Python list literals are interpreted as SymPy matrices.
 
     .. warning::
-        Note that this function uses ``eval`` (because it calls
-        ``sympy.sympify``), and thus shouldn't be used on unsanitized input.
+        This function uses ``eval`` (because it calls ``sympy.sympify``), and
+        thus should not be used on unsanitized input.
+
+    If `expr` is already a SymPy expression or a SymPy matrix, both the keys
+    and the values of `subs` are sympified (with the above special rules in
+    force) and used to to perform substitution using the ``.subs`` method of
+    `expr`.
 
-    If `expr` is a SymPy expression or a SymPy matrix, both the keys and values
-    of `subs` are sympified (without the above special rules) and used to to
-    perform substitution using the ``.subs`` method of `expr`.  The input
-    expression is not altered otherwise.
+    .. note::
+        Any (part of) argument to this function that gets sympified (for
+        example the values of `subs`) may be already a SymPy object.  In this
+        case it is left unchanged.  In particular, the commutativity of its
+        terms is not altered.  This possibly confusing effect is demonstrated
+        in the last example below.
 
     Parameters
     ----------
     expr : str or SymPy expression
-        The expression to be converted to a SymPy object.
+        Expression to be converted to a SymPy object.
     subs : dict or ``None`` (default)
-        The namespace of substitutions to be performed on the input `expr`.
+        Namespace of substitutions to be performed on the input `expr`.
         If `expr` is a string, the keys must be strings as well.  Otherwise
         they may be any sympifyable object. The values must be sympifyable
         objects.
@@ -121,22 +125,23 @@ def sympify(expr, subs=None):
 
     Examples
     --------
-        >>> from kwant.continuum import sympify
         >>> sympify('k_x * A(x) * k_x + V(x)')
         k_x*A(x)*k_x + V(x)     # as valid sympy object
 
-        or
-
-        >>> from kwant.continuum import sympify
         >>> sympify('k_x**2 + V', subs={'V': 'V_0 + V(x)'})
         k_x**2 + V(x) + V_0
 
-        >>> from kwant.continuum import sympify
-        >>> sympify('k_x**2 * sigma_plus',
-        ...         subs={'sigma_plus': [[0, 2], [0, 0]]})
+        >>> subs = {'sigma_plus': [[0, 2], [0, 0]]}
+        >>> sympify('k_x**2 * sigma_plus', subs)
         Matrix([
         [0, 2*k_x**2],
         [0,        0]])
+
+        >>> sympify('k_x * A(c) * k_x', subs={'c': 'x'})
+        k_x*A(x)*k_x
+        >>> sympify('k_x * A(c) * k_x', subs={'c': sympy.Symbol('x')})
+        A(x)*k_x**2
+
     """
     stored_value = None
     sympified_types = (sympy.Expr, sympy.matrices.MatrixBase)
diff --git a/kwant/continuum/discretizer.py b/kwant/continuum/discretizer.py
index c80c370f38a623dd3976548d2e334c54b285ecc6..bc375dc68a9c6ed651d6f9995530608044162298 100644
--- a/kwant/continuum/discretizer.py
+++ b/kwant/continuum/discretizer.py
@@ -88,34 +88,36 @@ def discretize(hamiltonian, coords=None, *, grid_spacing=1,
     `~kwant.continuum.discretize_symbolic` and feeding its result into
     `~kwant.continuum.build_discretized`.
 
+    .. warning::
+        This function uses ``eval`` (because it calls ``sympy.sympify``), and
+        thus should not be used on unsanitized input.
+
     Parameters
     ----------
-    hamiltonian : sympy.Expr or sympy.Matrix, or string
-        Symbolic representation of a continous Hamiltonian. When providing
-        a sympy expression, it is recommended to use ``momentum_operators``
-        defined in `~kwant.continuum` in order to provide
-        proper commutation properties. If a string is provided it will be
-        converted to a sympy expression using `kwant.continuum.sympify`.
+    hamiltonian : string, or sympy.Expr, or sympy.Matrix
+        Symbolic representation of a continous Hamiltonian.  It is
+        converted to a SymPy expression using `kwant.continuum.sympify`.
     coords : sequence of strings, or ``None`` (default)
         The coordinates for which momentum operators will be treated as
         differential operators. May contain only "x", "y" and "z" and must be
         sorted.  If not provided, `coords` will be obtained from the input
         Hamiltonian by reading the present coordinates and momentum operators.
     grid_spacing : int or float, default: 1
-        Grid spacing for the template Builder.
-    subs : dict, defaults to empty
-        A namespace of substitutions to be performed on the input
-        ``hamiltonian``. Values can be either strings or ``sympy`` objects.
-        It can be used to simplify input of matrices or alternate input before
-        proceeding further. For example:
+        Spacing of the (quadratic or cubic) discretization grid.
+    subs : dict or ``None`` (default)
+        A namespace of substitutions to be performed on `hamiltonian`. See
+        `kwant.continuum.sympify` for details. May be used to simplify input of
+        matrices or alternate input before proceeding further. For example:
         ``subs={'k': 'k_x + I * k_y'}`` or
-        ``subs={'s_z': [[1, 0], [0, -1]]}``.
+        ``subs={'sigma_plus': [[0, 2], [0, 0]]}``.
 
     Returns
     -------
-    template: `~kwant.builder.Builder`
+    model: `~kwant.builder.Builder`
         The translationally symmetric builder that corresponds to the provided
-        Hamiltonian.
+        Hamiltonian.  It may be passed directly into the method
+        `~kwant.builder.Builder.fill` of a builder to construct a system of a
+        desired shape.
     """
     tb, coords = discretize_symbolic(hamiltonian, coords, subs=subs)
 
@@ -128,30 +130,30 @@ def discretize_symbolic(hamiltonian, coords=None, *, subs=None):
     The two objects returned by this function may be used directly as the first
     two arguments for `~kwant.continuum.build_discretized`.
 
+    .. warning::
+        This function uses ``eval`` (because it calls ``sympy.sympify``), and
+        thus should not be used on unsanitized input.
+
     Parameters
     ----------
-    hamiltonian : sympy.Expr or sympy.Matrix, or string
-        Symbolic representation of a continuous Hamiltonian. When providing
-        a sympy expression, it is recommended to use ``momentum_operators``
-        defined in `~kwant.continuum` in order to provide
-        proper commutation properties. If a string is provided it will be
+    hamiltonian : string, or sympy.Expr, or sympy.Matrix
+        Symbolic representation of a continous Hamiltonian.  It is
         converted to a sympy expression using `kwant.continuum.sympify`.
     coords : sequence of strings, or ``None`` (default)
         The coordinates for which momentum operators will be treated as
         differential operators. May contain only "x", "y" and "z" and must be
         sorted.  If not provided, `coords` will be obtained from the input
         Hamiltonian by reading the present coordinates and momentum operators.
-    subs: dict, defaults to empty
-        A namespace of substitutions to be performed on the input
-        ``hamiltonian``. Values can be either strings or ``sympy`` objects.
-        It can be used to simplify input of matrices or alternate input before
-        proceeding further. For example:
+    subs : dict, defaults to empty
+        A namespace of substitutions to be performed on `hamiltonian`. See
+        `kwant.continuum.sympify` for details. May be used to simplify input of
+        matrices or alternate input before proceeding further. For example:
         ``subs={'k': 'k_x + I * k_y'}`` or
-        ``subs={'s_z': [[1, 0], [0, -1]]}``.
+        ``subs={'sigma_plus': [[0, 2], [0, 0]]}``.
 
     Returns
     -------
-    tb_hamiltonian: dict
+    tb_hamiltonian : dict
         Keys are tuples of integers; the offsets of the hoppings ((0, 0, 0) for
         the onsite). Values are symbolic expressions for the hoppings/onsite.
     coords : list of strings
@@ -213,39 +215,44 @@ def discretize_symbolic(hamiltonian, coords=None, *, subs=None):
     return tb, coords
 
 
-def build_discretized(tb_hamiltonian, coords, *,
-                      grid_spacing=1, subs=None):
-    """Create a template Builder from a symbolic tight-binding Hamiltonian.
+def build_discretized(tb_hamiltonian, coords, *, grid_spacing=1, subs=None):
+    """Create a template builder from a symbolic tight-binding Hamiltonian.
 
-    This return values of `~kwant.continuum.discretize_symbolic` may be used
+    The return values of `~kwant.continuum.discretize_symbolic` may be used
     directly for the first two arguments of this function.
 
+    .. warning::
+        This function uses ``eval`` (because it calls ``sympy.sympify``), and
+        thus should not be used on unsanitized input.
+
     Parameters
     ----------
     tb_hamiltonian : dict
-        Keys are tuples of integers: the offsets of the hoppings
-        ((0, 0, 0) for the onsite). Values are symbolic expressions
-        for the hoppings/onsite or expressions that can by sympified using
+        Keys must be the offsets of the hoppings, represented by tuples of
+        integers ((0, 0, 0) for onsite). Values must be symbolic expressions
+        for the hoppings/onsite or expressions that can by sympified with
         `kwant.continuum.sympify`.
     coords : sequence of strings
         The coordinates for which momentum operators will be treated as
         differential operators. May contain only "x", "y" and "z" and must be
         sorted.
     grid_spacing : int or float, default: 1
-        Grid spacing for the template Builder.
-    subs: dict, defaults to empty
-        A namespace of substitutions to be performed on the values of input
-        `tb_hamiltonian`. Values can be either strings or ``sympy`` objects.
-        It can be used to simplify input of matrices or alternate input before
-        proceeding further. For example:
-        ``subs={'k': 'k_x + I * k_y'}`` or
-        ``subs={'s_z': [[1, 0], [0, -1]]}``.
+        Spacing of the (quadratic or cubic) discretization grid.
+    subs : dict, defaults to empty
+        A namespace of substitutions to be performed on the values of
+        `tb_hamiltonian`. See `kwant.continuum.sympify` for details. May be
+        used to simplify input of matrices or alternate input before proceeding
+        further. For example: ``subs={'k': 'k_x + I * k_y'}`` or
+        ``subs={'sigma_plus': [[0, 2], [0, 0]]}``.
 
     Returns
     -------
-    template: `~kwant.builder.Builder`
+    model : `~kwant.builder.Builder`
         The translationally symmetric builder that corresponds to the provided
-        Hamiltonian.
+        Hamiltonian.  It may be passed directly into the method
+        `~kwant.builder.Builder.fill` of a builder to construct a system of a
+        desired shape.
+
     """
     if len(coords) == 0:
         raise ValueError('Discrete coordinates cannot be empty.')