diff --git a/brian2/stateupdaters/exponential_euler.py b/brian2/stateupdaters/exponential_euler.py
index 8dd9a1f83..7e967618d 100644
--- a/brian2/stateupdaters/exponential_euler.py
+++ b/brian2/stateupdaters/exponential_euler.py
@@ -1,9 +1,10 @@
 import sympy as sp
 
-from brian2.stateupdaters.base import StateUpdateMethod
+from .base import StateUpdateMethod
 
 __all__ = ['exponential_euler']
 
+
 def get_conditionally_linear_system(eqs):
     '''
     Convert equations into a linear system using sympy.
@@ -15,12 +16,28 @@ def get_conditionally_linear_system(eqs):
     
     Returns
     -------
-    
+    coefficients : dict of (sympy expression, sympy expression) tuples
+        For every variable x, a tuple (M, B) containing the coefficients M and
+        B (as sympy expressions) for M * x + B
     
     Raises
     ------
     ValueError
-        If the equations cannot be converted into an M * X + B form.
+        If one of the equations cannot be converted into a M * x + B form.
+
+    Examples
+    --------
+    >>> from brian2 import Equations
+    >>> eqs = Equations("""
+    ... dv/dt = (-v + w**2) / tau : 1
+    ... dw/dt = -w / tau : 1
+    ... """)
+    >>> system = get_conditionally_linear_system(eqs)
+    >>> print(system['v'])
+    (-1/tau, w**2/tau)
+    >>> print(system['w'])
+    (-1/tau, 0)
+
     '''
     diff_eqs = eqs.substituted_expressions
     
@@ -48,8 +65,22 @@ def get_conditionally_linear_system(eqs):
     
     return coefficients
 
+
 class ExponentialEulerStateUpdater(StateUpdateMethod):
+    '''
+    A state updater for conditionally linear equations, i.e. equations where
+    each variable only depends linearly on itself (but possibly non-linearly
+    on other variables). Typical Hodgkin-Huxley equations fall into this
+    category, it is therefore the default integration method used in the
+    GENESIS simulator, for example.
+    '''
     def can_integrate(self, equations, namespace, specifiers):
+        '''
+        Return whether the given equations can be integrated using this
+        state updater. This method tests for conditional linearity by
+        executing `get_conditionally_linear_system` and returns ``True`` in
+        case this call succeeds.
+        '''
         if equations.is_stochastic:
             return False
         
@@ -87,5 +118,7 @@ def __call__(self, equations, namespace=None, specifiers=None):
             code += ['{var} = _{var}'.format(var=var)]
             
         return '\n'.join(code)
+    # Copy doc from parent class
+    __call__.__doc__ = StateUpdateMethod.__call__.__doc__
 
 exponential_euler = ExponentialEulerStateUpdater() 
\ No newline at end of file
diff --git a/examples/I-F_curve_Hodgkin_Huxley.py b/examples/I-F_curve_Hodgkin_Huxley.py
new file mode 100644
index 000000000..2a41c262c
--- /dev/null
+++ b/examples/I-F_curve_Hodgkin_Huxley.py
@@ -0,0 +1,52 @@
+'''
+Input-Frequency curve of a HH model
+Network: 100 unconnected Hodgin-Huxley neurons with an input current I.
+The input is set differently for each neuron.
+
+This simulation should use exponential Euler integration
+'''
+from pylab import *
+from brian2 import *
+
+BrianLogger.log_level_info()
+
+language = PythonLanguage()
+#language = CPPLanguage()
+
+N = 100
+
+# Parameters
+area = 20000 * umetre ** 2
+Cm = (1 * ufarad * cm ** -2) * area
+gl = (5e-5 * siemens * cm ** -2) * area
+El = -65 * mV
+EK = -90 * mV
+ENa = 50 * mV
+g_na = (100 * msiemens * cm ** -2) * area
+g_kd = (30 * msiemens * cm ** -2) * area
+VT = -63 * mV
+
+# The model
+eqs = Equations('''
+dv/dt = (gl*(El-v) - g_na*(m*m*m)*h*(v-ENa) - g_kd*(n*n*n*n)*(v-EK) + I)/Cm : volt
+dm/dt = 0.32*(mV**-1)*(13.*mV-v+VT)/
+    (exp((13.*mV-v+VT)/(4.*mV))-1.)/ms*(1-m)-0.28*(mV**-1)*(v-VT-40.*mV)/
+    (exp((v-VT-40.*mV)/(5.*mV))-1.)/ms*m : 1
+dn/dt = 0.032*(mV**-1)*(15.*mV-v+VT)/
+    (exp((15.*mV-v+VT)/(5.*mV))-1.)/ms*(1.-n)-.5*exp((10.*mV-v+VT)/(40.*mV))/ms*n : 1
+dh/dt = 0.128*exp((17.*mV-v+VT)/(18.*mV))/ms*(1.-h)-4./(1+exp((40.*mV-v+VT)/(5.*mV)))/ms*h : 1
+I : amp
+''')
+# Threshold and refractoriness are only used for spike counting
+group = NeuronGroup(N, equations=eqs, threshold='is_active * (v > -40*mV)',
+                    language=language)
+group.refractory = 1*ms
+group.v = El
+group.I = linspace(0 * nA, 0.7 * nA, N)
+
+monitor = SpikeMonitor(group)
+
+duration = 2 * second
+run(duration)
+plot(group.I / nA, monitor.count / duration)
+show()
diff --git a/examples/I-F_curve.py b/examples/I-F_curve_LIF.py
similarity index 96%
rename from examples/I-F_curve.py
rename to examples/I-F_curve_LIF.py
index 5435307c4..c932a7103 100644
--- a/examples/I-F_curve.py
+++ b/examples/I-F_curve_LIF.py
@@ -22,7 +22,7 @@
 
 monitor = SpikeMonitor(group)
 
-duration = 5 * second
+duration = 1 * second
 run(duration)
 plot(group.v0 / mV, monitor.count / duration)
 show()