diff --git a/CHANGES.rst b/CHANGES.rst
index d9fb7886..e580d340 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -1,3 +1,10 @@
+v0.4.1
+======
+- Fixes for enhanced robustness:
+  - ``.kinetics.ode.get_odesys``
+  - ``.chemistry.as_per_substance_array``
+- Minor changes.
+
 v0.4.0
 ======
 - Multiple fixes throughout
diff --git a/chempy/chemistry.py b/chempy/chemistry.py
index 30e1942a..43718e78 100644
--- a/chempy/chemistry.py
+++ b/chempy/chemistry.py
@@ -1301,7 +1301,7 @@ def as_per_substance_array(self, cont, dtype='float64', unit=None, raise_on_unk=
                         raise KeyError("Unkown substance key: %s" % k)
             cont = [cont[k] for k in substance_keys]
 
-        cont = np.asarray(cont, dtype=dtype)
+        cont = np.atleast_1d(np.asarray(cont, dtype=dtype).squeeze())
         if cont.shape[-1] != self.ns:
             raise ValueError("Incorrect size")
         return cont*(unit if unit is not None else 1)
diff --git a/chempy/kinetics/_rates.py b/chempy/kinetics/_rates.py
index ca79b559..b8b10be5 100644
--- a/chempy/kinetics/_rates.py
+++ b/chempy/kinetics/_rates.py
@@ -16,7 +16,7 @@ class TPolyMassAction(TPoly, MassAction):
     argument_names = ('temperature_offset', Ellipsis)
     parameter_keys = TPoly.parameter_keys
 
-    def rate_coeff(self, variables, backend):
+    def rate_coeff(self, variables, backend=math):
         return self.eval_poly(variables, backend)
 
 
@@ -27,7 +27,7 @@ def rate_coeff(self, variables, backend):
 class PiecewiseTPolyMassAction(TPiecewisePoly, MassAction):
     parameter_keys = TPoly.parameter_keys
 
-    def rate_coeff(self, variables, backend):
+    def rate_coeff(self, variables, backend=math):
         return self.eval_poly(variables, backend)
 
 
@@ -39,19 +39,27 @@ def g_value(self, variables, backend):
         return self.eval_poly(variables, backend)
 
 
+class RTPolyRadiolytic(RTPoly, Radiolytic):
+    nargs = None
+    parameter_keys = Radiolytic.parameter_keys + RTPoly.parameter_keys
+
+    def g_value(self, variables, backend):
+        return self.eval_poly(variables, backend)
+
+
 class RTPolyMassAction(RTPoly, MassAction):
     """ Arguments: temperature_offset, c0, c1, ... """
     parameter_keys = RTPoly.parameter_keys
     nargs = None
 
-    def rate_coeff(self, variables, backend):
+    def rate_coeff(self, variables, backend=math):
         return self.eval_poly(variables, backend)
 
 
 class _Log10XPolyMassAction(MassAction):
     skip_poly = 1  # kunit
 
-    def rate_coeff(self, variables, backend):
+    def rate_coeff(self, variables, backend=math):
         k_unit = self.arg(variables, 0)
         return 10**self.eval_poly(variables, backend)*k_unit
 
diff --git a/chempy/kinetics/ode.py b/chempy/kinetics/ode.py
index d1a8d218..894c222f 100644
--- a/chempy/kinetics/ode.py
+++ b/chempy/kinetics/ode.py
@@ -162,7 +162,6 @@ def get_odesys(rsys, include_params=True, substitutions=None,
             if ratex.unique_keys is not None:
                 unique_keys.extend(ratex.unique_keys)
                 p_defaults.extend(ratex.args)
-
     if unit_registry is None:
         def pre_processor(x, y, p):
             return (
@@ -184,9 +183,9 @@ def post_processor(x, y, p):
         p_units = [get_derived_unit(unit_registry, k) for k in param_keys]
         new_r_exprs = []
         for ratex in r_exprs:
-            _pu, _new_rates = ratex.dedimensionalisation(unit_registry)
+            _pu, _new_rate = ratex._recursive_as_RateExpr().dedimensionalisation(unit_registry)
             p_units.extend(_pu)
-            new_r_exprs.append(_new_rates)
+            new_r_exprs.append(_new_rate)
         r_exprs = new_r_exprs
 
         time_unit = get_derived_unit(unit_registry, 'time')
diff --git a/chempy/kinetics/rates.py b/chempy/kinetics/rates.py
index bd8d56f3..21e6bcda 100644
--- a/chempy/kinetics/rates.py
+++ b/chempy/kinetics/rates.py
@@ -31,6 +31,22 @@ def rxn(self, value):
             if isinstance(arg, RateExpr):
                 arg.rxn = value
 
+    def _recursive_as_RateExpr(self):
+        new_args = []
+        for arg in self.args:
+            if isinstance(arg, Expr):
+                new_args.append(arg)
+            else:
+                if hasattr(arg, '_as_RateExpr'):
+                    new_args.append(arg._as_RateExpr(self.rxn))
+                else:
+                    new_args.append(arg)
+        if self.kw is None:
+            kw = {}
+        else:
+            kw = {k: getattr(self, k) for k in self.kw}
+        return self.__class__(new_args, self.unique_keys, **kw)
+
     @classmethod
     @deprecated(use_instead=Expr.from_callback)
     def subclass_from_callback(cls, cb, cls_attrs=None):
@@ -91,7 +107,7 @@ class _Radiolytic(RadiolyticBase):
         parameter_keys = (doserate_name, 'density')
         print_name = 'Radiolytic' if doserate_name == 'doserate' else ('Radiolytic{'+doserate_name+'}')
 
-        def g_value(self, variables, backend):  # for subclasses
+        def g_value(self, variables, backend=math):  # for subclasses
             return self.arg(variables, 0, backend=backend)
 
         def __call__(self, variables, backend=math):
@@ -168,6 +184,22 @@ def as_mass_action(self, variables, backend=math):
 
 
 class ArrheniusMassAction(MassAction):
+    """ Rate expression for a Arrhenius-type of rate
+
+    Examples
+    --------
+    >>> from math import exp
+    >>> from chempy import Reaction
+    >>> from chempy.units import allclose, default_units as u
+    >>> A = 1e11 / u.second
+    >>> Ea_over_R = 42e3/8.3145 * u.K**-1
+    >>> ratex = ArrheniusMassAction([A, Ea_over_R])
+    >>> rxn = Reaction({'R'}, {'P'}, ratex)
+    >>> dRdt = rxn.rate({'R': 3*u.M, 'temperature': 298.15*u.K})['R']
+    >>> allclose(dRdt, -3*1e11*exp(-42e3/8.3145/298.15)*u.M/u.s)
+    True
+
+    """
     argument_names = ('A', 'Ea_over_R')
     parameter_keys = ('temperature',)
 
diff --git a/chempy/tests/test_chemistry.py b/chempy/tests/test_chemistry.py
index cb7d5366..95ef28c5 100644
--- a/chempy/tests/test_chemistry.py
+++ b/chempy/tests/test_chemistry.py
@@ -187,8 +187,7 @@ def test_ReactionSystem__as_per_substance_array_dict():
     m = default_units.metre
     M = default_units.molar
     rs = ReactionSystem([], [Substance('H2O')])
-    c = rs.as_per_substance_array({'H2O': 1*M},
-                                  unit=M)
+    c = rs.as_per_substance_array({'H2O': 1*M}, unit=M)
     assert c.dimensionality == M.dimensionality
     assert abs(c[0]/(1000*mol/m**3) - 1) < 1e-16
 
diff --git a/chempy/tests/test_units.py b/chempy/tests/test_units.py
index 41bb4699..74efd242 100644
--- a/chempy/tests/test_units.py
+++ b/chempy/tests/test_units.py
@@ -11,9 +11,9 @@
 
 from ..util.testing import requires
 from ..units import (
-    allclose, get_derived_unit, is_unitless, linspace,
+    allclose, get_derived_unit, is_unitless, linspace, logspace_from_lin,
     SI_base_registry, unitless_in_registry, format_string, get_physical_quantity,
-    to_unitless, magnitude, default_unit_in_registry, Backend,
+    to_unitless, magnitude, default_unit_in_registry, Backend, latex_of_unit,
     unit_of, unit_registry_to_human_readable, units_library,
     unit_registry_from_human_readable, _sum, UncertainQuantity,
     default_units as u
@@ -104,6 +104,7 @@ def test_to_unitless():
     assert abs(to_unitless(3/(u.second*u.molar),
                            u.metre**3/u.mole/u.second) - 3e-3) < 1e-12
     assert abs(to_unitless(2*u.dm3, u.cm3) - 2000) < 1e-12
+    assert abs(to_unitless(2*u.m3, u.dm3) - 2000) < 1e-12
     assert (float(to_unitless(UncertainQuantity(2, u.dm3, .3), u.cm3)) - 2000) < 1e-12
 
     g1 = UncertainQuantity(4.46, u.per100eV, 0)
@@ -137,6 +138,13 @@ def test_linspace():
     assert abs(to_unitless(ls[0], u.hour) - 2/3600.) < 1e-15
 
 
+@requires(units_library)
+def test_logspace_from_lin():
+    ls = logspace_from_lin(2*u.second, 3*u.second)
+    assert abs(to_unitless(ls[0], u.hour) - 2/3600.) < 1e-15
+    assert abs(to_unitless(ls[-1], u.hour) - 3/3600.) < 1e-15
+
+
 @requires(units_library)
 def test_get_derived_unit():
     registry = SI_base_registry.copy()
@@ -323,3 +331,8 @@ def test_joule_html():
     joule_htm = 'kg&sdot;m<sup>2</sup>/s<sup>2</sup>'
     joule = u.J.dimensionality.simplified
     assert joule.html == joule_htm
+
+
+@requires(units_library)
+def test_latex_of_unit():
+    assert latex_of_unit(u.gram/u.metre**2) == r'\mathrm{\frac{g}{m^{2}}}'
diff --git a/chempy/units.py b/chempy/units.py
index d47888e3..1dde2d2a 100644
--- a/chempy/units.py
+++ b/chempy/units.py
@@ -5,19 +5,21 @@
 - ``chempy.units.default_constants``
 - ``chempy.units.SI_base_registry``
 
+together with some functions.
+
 Currently `quantities <https://pypi.python.org/pypi/quantities>`_ is used as
 the underlying package to handle units. If it is possible you should try to
-only use the `chempy.units` module (in case ``ChemPy`` changes this backend).
+only use the `chempy.units` module (in case ``ChemPy`` changes this backend),
+and avoid relying on any attributes of the Quantity instances (and rather use
+functions in `chempy.units`).
 
 """
 from __future__ import (absolute_import, division, print_function)
 
-from operator import mul
 from functools import reduce
+from operator import mul
 
 from ._util import NameSpace
-# Currently we use quantities for units. This may change, therefore use this
-# file for all units. A requirement is first-class numpy support.
 
 units_library = 'quantities'  # info used for selective testing.
 
@@ -40,6 +42,7 @@
     default_units = NameSpace(pq)
     default_units.decimetre = pq.UnitQuantity(
         'decimetre',  default_units.m / 10.0, u_symbol='dm')
+    default_units.m3 = default_units.metre**3
     default_units.dm3 = default_units.decimetre**3
     default_units.cm3 = default_units.centimetre**3
     if not hasattr(default_units, 'molar'):
@@ -147,6 +150,18 @@ def unit_registry_to_human_readable(unit_registry):
     return new_registry
 
 
+def latex_of_unit(quant):
+    """ Returns LaTeX reperesentation of the unit of a quantity
+
+    Examples
+    --------
+    >>> print(latex_of_unit(1/default_units.kelvin))
+    \\mathrm{\\frac{1}{K}}
+
+    """
+    return quant.dimensionality.latex.strip('$')
+
+
 def unit_registry_from_human_readable(unit_registry):
     """ Deserialization of unit_registry. """
     if unit_registry is None:
@@ -308,7 +323,15 @@ def allclose(a, b, rtol=1e-8, atol=None):
 
 
 def linspace(start, stop, num=50):
-    """ Analogous to ``numpy.linspace``. """
+    """ Analogous to ``numpy.linspace``.
+
+    Examples
+    --------
+    >>> abs(linspace(2, 8, num=3)[1] - 5) < 1e-15
+    True
+
+    """
+
     # work around for quantities v0.10.1 and NumPy
     import numpy as np
     unit = unit_of(start)
@@ -317,6 +340,22 @@ def linspace(start, stop, num=50):
     return np.linspace(start_, stop_, num)*unit
 
 
+def logspace_from_lin(start, stop, num=50):
+    """ Logarithmically spaced data points
+
+    Examples
+    --------
+    >>> abs(logspace_from_lin(2, 8, num=3)[1] - 4) < 1e-15
+    True
+
+    """
+    import numpy as np
+    unit = unit_of(start)
+    start_ = np.log2(to_unitless(start, unit))
+    stop_ = np.log2(to_unitless(stop, unit))
+    return np.exp2(np.linspace(start_, stop_, num))*unit
+
+
 def _sum(iterable):
     try:
         result = next(iterable)
diff --git a/chempy/util/_expr.py b/chempy/util/_expr.py
index d8bebebf..647df07b 100644
--- a/chempy/util/_expr.py
+++ b/chempy/util/_expr.py
@@ -12,6 +12,7 @@
 from functools import reduce
 from itertools import chain
 from operator import add, mul, truediv, sub
+from .pyutil import defaultkeydict
 
 
 class Expr(object):
@@ -55,8 +56,8 @@ class Expr(object):
     >>> cv = {s.name: EinsteinSolid([einT(s), s.mass]) for s in (Al, Be)}
     >>> print('%.4f' % cv['Al']({'temperature': 273.15, 'molar_gas_constant': 8.3145}))  # J/(g*K)
     0.8108
-    >>> import sympy
-    >>> print(cv['Be']({'temperature': sympy.Symbol('T'), 'molar_gas_constant': sympy.Symbol('R')}, backend=sympy))
+    >>> import sympy; from sympy import Symbol as Symb
+    >>> print(cv['Be']({'temperature': Symb('T'), 'molar_gas_constant': Symb('R')}, backend=sympy))
     112105.346283965*R/(T**2*sinh(580.32/T)**2)
 
     Attributes
@@ -104,6 +105,44 @@ def __init__(self, args, unique_keys=None, **kwargs):
         if kwargs:
             raise ValueError("Unexpected keyword arguments %s" % kwargs)
 
+    @classmethod
+    def from_callback(cls, callback, **kwargs):
+        """ Factory of subclasses
+
+        Parameters
+        ----------
+        callback : callable
+            signature: *args, backend=None
+        argument_names : tuple of str, optional
+        parameter_keys : tuple of str, optional,
+        kw : dict, optional
+        nargs : int, optional
+
+        Examples
+        --------
+        >>> from operator import add; from functools import reduce
+        >>> def poly(args, x, backend=None):
+        ...     x0 = args[0]
+        ...     return reduce(add, [c*(x-x0)**i for i, c in enumerate(args[1:])])
+        ...
+        >>> Poly = Expr.from_callback(poly, parameter_keys=('x',), argument_names=('x0', Ellipsis))
+        >>> p = Poly([1, 3, 2, 5])
+        >>> p({'x': 7}) == 3 + 2*(7-1) + 5*(7-1)**2
+        True
+        >>> q = Poly([1, 3, 2, 5], unique_keys=('x0_q',))
+        >>> q({'x': 7, 'x0_q': 0}) == 3 + 2*7 + 5*7**2
+        True
+
+        """
+        class Wrapper(cls):
+            def __call__(self, variables, backend=math):
+                args = self.all_args(variables, backend=backend)
+                params = self.all_params(variables, backend=backend)
+                return callback(args, *params, backend=backend)
+        for k, v in kwargs.items():
+            setattr(Wrapper, k, v)
+        return Wrapper
+
     def __call__(self, variables, backend=None):
         raise NotImplementedError("Subclass and implement __call__")
 
@@ -134,20 +173,20 @@ def __repr__(self):
     def string(self, arg_fmt=str, **kwargs):
         return self._str(arg_fmt, **kwargs)
 
-    def arg(self, variables, index, backend=None):
+    def arg(self, variables, index, backend=None, evaluate=True):
         if isinstance(index, str):
             index = self.argument_names.index(index)
         if self.unique_keys is None or len(self.unique_keys) <= index:
             res = self.args[index]
         else:
             res = variables.get(self.unique_keys[index], self.args[index])
-        if isinstance(res, Expr):
+        if isinstance(res, Expr) and evaluate:
             return res(variables, backend=backend)
         else:
             return res
 
-    def all_args(self, variables, backend=None):
-        return [self.arg(variables, i, backend) for i in range(len(self.args))]
+    def all_args(self, variables, backend=None, evaluate=True):
+        return [self.arg(variables, i, backend, evaluate) for i in range(len(self.args))]
 
     def all_params(self, variables, backend=None):
         return [v(variables, backend=backend) if isinstance(v, Expr) else v for v
@@ -189,9 +228,9 @@ def dedimensionalisation(self, unit_registry, variables={}, backend=math):
         """
         from ..units import default_unit_in_registry, to_unitless
         units = [None if isinstance(arg, Expr) else default_unit_in_registry(arg, unit_registry) for arg
-                 in self.all_args(variables, backend=backend)]
+                 in self.all_args(variables, backend=backend, evaluate=False)]
         new_units, unitless_args = [], []
-        for arg, unit in zip(self.all_args(variables, backend=backend), units):
+        for arg, unit in zip(self.all_args(variables, backend=backend, evaluate=False), units):
             if isinstance(arg, Expr):
                 if unit is not None:
                     raise ValueError()
@@ -206,12 +245,14 @@ def dedimensionalisation(self, unit_registry, variables={}, backend=math):
             kw = {k: getattr(self, k) for k in self.kw}
         return new_units, self.__class__(unitless_args, self.unique_keys, **kw)
 
-    def _sympy_format(self, method, variables, backend):
+    def _sympy_format(self, method, variables, backend, default):
         variables = variables or {}
         if backend is None:
             import sympy as backend
-        variables = {k: v if isinstance(v, Expr) else backend.Symbol(v) for k, v in variables.items()}
-        expr = self(variables, backend=backend)
+        variables = defaultkeydict(
+            None if default is None else (lambda k: backend.Symbol(default(k))),
+            {k: v if isinstance(v, Expr) else backend.Symbol(v) for k, v in variables.items()})
+        expr = self(variables, backend=backend).simplify()
         if method == 'latex':
             return backend.latex(expr)
         elif method == 'unicode':
@@ -222,12 +263,14 @@ def _sympy_format(self, method, variables, backend):
         else:
             raise NotImplementedError("Unknown method: %s" % method)
 
-    def latex(self, variables=None, backend=None):
+    def latex(self, variables=None, backend=None, default=None):
         r"""
         Parameters
         ----------
         variables : dict
         backend : module
+        default : callable
+            Format string based on missing key, signature: str -> str.
 
         Examples
         --------
@@ -243,45 +286,7 @@ def latex(self, variables=None, backend=None):
         Requires SymPy
 
         """
-        return self._sympy_format('latex', variables, backend)
-
-    @classmethod
-    def from_callback(cls, callback, **kwargs):
-        """ Factory of subclasses
-
-        Parameters
-        ----------
-        callback : callable
-            signature: *args, backend=None
-        argument_names : tuple of str, optional
-        parameter_keys : tuple of str, optional,
-        kw : dict, optional
-        nargs : int, optional
-
-        Examples
-        --------
-        >>> from operator import add; from functools import reduce
-        >>> def poly(args, x, backend=None):
-        ...     x0 = args[0]
-        ...     return reduce(add, [c*(x-x0)**i for i, c in enumerate(args[1:])])
-        ...
-        >>> Poly = Expr.from_callback(poly, parameter_keys=('x',), argument_names=('x0', Ellipsis))
-        >>> p = Poly([1, 3, 2, 5])
-        >>> p({'x': 7}) == 3 + 2*(7-1) + 5*(7-1)**2
-        True
-        >>> q = Poly([1, 3, 2, 5], unique_keys=('x0_q',))
-        >>> q({'x': 7, 'x0_q': 0}) == 3 + 2*7 + 5*7**2
-        True
-
-        """
-        class Wrapper(cls):
-            def __call__(self, variables, backend=math):
-                args = self.all_args(variables, backend=backend)
-                params = self.all_params(variables, backend=backend)
-                return callback(args, *params, backend=backend)
-        for k, v in kwargs.items():
-            setattr(Wrapper, k, v)
-        return Wrapper
+        return self._sympy_format('latex', variables, backend=backend, default=default)
 
     def __eq__(self, other):
         if self.__class__ != other.__class__:
diff --git a/chempy/util/bkh.py b/chempy/util/bkh.py
index 4365f496..767c9296 100644
--- a/chempy/util/bkh.py
+++ b/chempy/util/bkh.py
@@ -8,13 +8,13 @@
 from itertools import chain
 
 from chempy.kinetics.ode import get_odesys
-from chempy.units import to_unitless, linspace
+from chempy.units import to_unitless, linspace, logspace_from_lin
 
 
 def integration_with_sliders(
         rsys, tend, c0, parameters, fig_kwargs=None, unit_registry=None, output_conc_unit=None,
         output_time_unit=None, slider_kwargs=None, x_axis_type="linear", y_axis_type="linear",
-        integrate_kwargs=None):
+        integrate_kwargs=None, substitutions=None):
 
     import numpy as np
     from bokeh.plotting import Figure
@@ -24,7 +24,7 @@ def integration_with_sliders(
     if slider_kwargs is None:
         slider_kwargs = {}
     odesys, state_keys, rarg_keys, p_units = get_odesys(
-        rsys, unit_registry=unit_registry,
+        rsys, unit_registry=unit_registry, substitutions=substitutions,
         output_conc_unit=output_conc_unit,
         output_time_unit=output_time_unit)[:4]
     if output_conc_unit is None:
@@ -33,7 +33,13 @@ def integration_with_sliders(
         output_conc_unit = 1
 
     param_keys = list(chain(state_keys, rarg_keys))
-    tout = linspace(tend*0, tend)
+    if x_axis_type == 'linear':
+        tout = linspace(tend*0, tend)
+    elif x_axis_type == 'log':
+        tout = logspace_from_lin(tend*1e-9, tend)
+    else:
+        raise NotImplementedError("Unknown x_axis_type: %s" % x_axis_type)
+
     tout, Cout, info = odesys.integrate(tout, c0, parameters, **(integrate_kwargs or {}))
     sources = [ColumnDataSource(data={
         'tout': to_unitless(tout, output_time_unit),
@@ -73,7 +79,7 @@ def _dict_to_unitless(d, u):
                        slider_kwargs.get(k, dict(start=v/10, end=v*10, step=v/10)),
                        u)))
         for k, v, u in zip(param_keys, p_ul, p_units)])
-    all_widgets = list(chain(c0_widgets.values(), param_widgets.values()))
+    all_widgets = list(chain(param_widgets.values(), c0_widgets.values()))
 
     def update_data(attrname, old, new):
         _c0 = defaultdict(lambda: 0*output_conc_unit)
diff --git a/chempy/util/pyutil.py b/chempy/util/pyutil.py
index 5d65ab10..7adf6052 100644
--- a/chempy/util/pyutil.py
+++ b/chempy/util/pyutil.py
@@ -4,7 +4,7 @@
 """
 from __future__ import (absolute_import, division, print_function)
 
-from collections import namedtuple, Mapping
+from collections import defaultdict, namedtuple, Mapping
 from functools import wraps
 import os
 import warnings
@@ -27,6 +27,29 @@ def wrapper(*args):
     return decorator
 
 
+class defaultkeydict(defaultdict):
+    """ defaultdict where default_factory should have the signature key -> value
+
+    Examples
+    --------
+    >>> d = defaultkeydict(lambda k: '[%s]' % k, {'a': '[a]', 'b': '[B]'})
+    >>> d['a']
+    '[a]'
+    >>> d['b']
+    '[B]'
+    >>> d['c']
+    '[c]'
+
+    """
+
+    def __missing__(self, key):
+        if self.default_factory is None:
+            raise KeyError("Missing key: %s" % key)
+        else:
+            self[key] = self.default_factory(key)
+        return self[key]
+
+
 def defaultnamedtuple(typename, field_names, defaults=()):
     """ Generates a new subclass of tuple with default values.
 
diff --git a/chempy/util/tests/test_expr.py b/chempy/util/tests/test_expr.py
index 1aeec8e6..5f320b6e 100644
--- a/chempy/util/tests/test_expr.py
+++ b/chempy/util/tests/test_expr.py
@@ -270,18 +270,25 @@ def test_Expr__latex():
     p = Poly([1, 2, 3, 4])
     import sympy
     t = sympy.Symbol('t')
-    ref = sympy.latex(2 + 3*(t-1) + 4*(t-1)**2)
+    ref = sympy.latex((2 + 3*(t-1) + 4*(t-1)**2).simplify())
     assert p.latex({'x': 't'}) == ref
 
     TE = Poly([3, 7, 5])
     cv_Al = _get_cv()['Al']
     T, E, R, m = sympy.symbols('T E R m')
     _TE = 7 + 5*(E-3)
-    ref = sympy.latex((3*R*(_TE/(2*T))**2 * sympy.sinh(_TE/(2*T))**-2)/m)
+    ref = sympy.latex(((3*R*(_TE/(2*T))**2 * sympy.sinh(_TE/(2*T))**-2)/m).simplify())
     cv_Al.unique_keys = ('TE_Al', 'm_Al')
     res = cv_Al.latex({'TE_Al': TE, 'temperature': 'T', 'x': 'E', 'molar_gas_constant': 'R', 'm_Al': 'm'})
     assert res == ref
 
+    X = sympy.symbols('X')
+    _TE2 = 7 + 5*(X-3)
+    ref2 = sympy.latex(((3*R*(_TE2/(2*T))**2 * sympy.sinh(_TE2/(2*T))**-2)/m).simplify())
+    res2 = cv_Al.latex({'TE_Al': TE, 'temperature': 'T', 'molar_gas_constant': 'R', 'm_Al': 'm'},
+                       default=lambda k: k.upper())
+    assert res2 == ref2
+
 
 def test_Expr__single_arg():
     class Pressure(Expr):
diff --git a/chempy/util/tests/test_pyutil.py b/chempy/util/tests/test_pyutil.py
index ebfae87f..2dfe8a21 100644
--- a/chempy/util/tests/test_pyutil.py
+++ b/chempy/util/tests/test_pyutil.py
@@ -1,7 +1,7 @@
 # -*- coding: utf-8 -*-
 from __future__ import (absolute_import, division, print_function)
 
-from ..pyutil import defaultnamedtuple
+from ..pyutil import defaultkeydict, defaultnamedtuple
 
 
 def test_defaultnamedtuple():
@@ -37,3 +37,8 @@ def __new__(cls, x2, y2=_default_y):
     p3 = MySubclass(9)
     assert p3.x == 3
     assert p3.y == 10**.5
+
+
+def test_defaultkeydict():
+    d = defaultkeydict(lambda k: k*2)
+    assert d['as'] == 'asas'
diff --git a/examples/bokeh_interactive_arrhenius_units.py b/examples/bokeh_interactive_arrhenius_units.py
index ebe865c8..94d8e002 100644
--- a/examples/bokeh_interactive_arrhenius_units.py
+++ b/examples/bokeh_interactive_arrhenius_units.py
@@ -18,7 +18,7 @@
     Af, Ab, Ea, Er = 1e16/u.molar/u.s, 1.5e15/u.s, 72e3*u.J/u.mol, -12e3*u.J/u.mol
     curdoc().add_root(integration_with_sliders(
         get_rsys(Af, Ab, Ea, Er), tend=3*u.s,
-        c0=defaultdict(float, {'Fe+3': 3e-3*u.molar, 'SCN-': 1.5e-3*u.molar}),
+        c0=defaultdict(lambda: 0*u.molar, {'Fe+3': 3e-3*u.molar, 'SCN-': 1.5e-3*u.molar}),
         parameters={'temperature': 298.15*u.K},
         slider_kwargs={'temperature': dict(start=273.15*u.K, end=313.15*u.K, step=.05*u.K)},
         unit_registry=SI_base_registry,