diff --git a/GPy/core/mapping.py b/GPy/core/mapping.py
index 8cfd9ad0..efd9476f 100644
--- a/GPy/core/mapping.py
+++ b/GPy/core/mapping.py
@@ -34,7 +34,7 @@ class Mapping(Parameterized):
         raise NotImplementedError
 
     def df_dtheta(self, dL_df, X):
-        """The gradient of the outputs of the mapping with respect to each of the parameters.
+        """The gradient of the outputs of the multi-layer perceptron with respect to each of the parameters.
 
         :param dL_df: gradient of the objective with respect to the function.
         :type dL_df: ndarray (num_data x output_dim)
@@ -50,7 +50,7 @@ class Mapping(Parameterized):
         """
         Plots the mapping associated with the model.
           - In one dimension, the function is plotted.
-          - In two dimensions, a contour-plot shows the function
+          - In two dimsensions, a contour-plot shows the function
           - In higher dimensions, we've not implemented this yet !TODO!
 
         Can plot only part of the data and part of the posterior functions
@@ -65,14 +65,6 @@ class Mapping(Parameterized):
         else:
             raise NameError, "matplotlib package has not been imported."
 
-class Bijective_mapping(Mapping):
-    """This is a mapping that is bijective, i.e. you can go from X to f and also back from f to X. The inverse mapping is called g()."""
-    def __init__(self, input_dim, output_dim, name='bijective_mapping'):
-        super(Bijective_apping, self).__init__(name=name)
-
-    def g(self, f):
-        """Inverse mapping from output domain of the function to the inputs."""
-        raise NotImplementedError
 
 from model import Model
 
diff --git a/GPy/core/parameterization/observable_array.py b/GPy/core/parameterization/observable_array.py
deleted file mode 100644
index fc9d6cf2..00000000
--- a/GPy/core/parameterization/observable_array.py
+++ /dev/null
@@ -1,137 +0,0 @@
-# Copyright (c) 2012, GPy authors (see AUTHORS.txt).
-# Licensed under the BSD 3-clause license (see LICENSE.txt)
-
-__updated__ = '2014-03-31'
-
-import numpy as np
-from parameter_core import Observable, Pickleable
-
-class ObsAr(np.ndarray, Pickleable, Observable):
-    """
-    An ndarray which reports changes to its observers.
-    The observers can add themselves with a callable, which
-    will be called every time this array changes. The callable
-    takes exactly one argument, which is this array itself.
-    """
-    __array_priority__ = -1 # Never give back ObsAr
-    def __new__(cls, input_array, *a, **kw):
-        if not isinstance(input_array, ObsAr):
-            obj = np.atleast_1d(np.require(input_array, dtype=np.float64, requirements=['W', 'C'])).view(cls)
-        else: obj = input_array
-        #cls.__name__ = "ObsAr" # because of fixed printing of `array` in np printing
-        super(ObsAr, obj).__init__(*a, **kw)
-        return obj
-
-    def __array_finalize__(self, obj):
-        # see InfoArray.__array_finalize__ for comments
-        if obj is None: return
-        self._observer_callables_ = getattr(obj, '_observer_callables_', None)
-
-    def __array_wrap__(self, out_arr, context=None):
-        return out_arr.view(np.ndarray)
-
-    def copy(self):
-        memo = {}
-        memo[id(self)] = self
-        return self.__deepcopy__(memo)
-
-    def __deepcopy__(self, memo):
-        s = self.__new__(self.__class__, input_array=self.view(np.ndarray).copy())
-        memo[id(self)] = s
-        import copy
-        s.__dict__.update(copy.deepcopy(self.__dict__, memo))
-        return s
-
-    def __reduce__(self):
-        func, args, state = super(ObsAr, self).__reduce__()
-        return func, args, (state, Pickleable.__getstate__(self))
-
-    def __setstate__(self, state):
-        np.ndarray.__setstate__(self, state[0])
-        Pickleable.__setstate__(self, state[1])
-
-    def __setitem__(self, s, val):
-        super(ObsAr, self).__setitem__(s, val)
-        self.notify_observers()
-
-    def __getslice__(self, start, stop):
-        return self.__getitem__(slice(start, stop))
-
-    def __setslice__(self, start, stop, val):
-        return self.__setitem__(slice(start, stop), val)
-
-    def __ilshift__(self, *args, **kwargs):
-        r = np.ndarray.__ilshift__(self, *args, **kwargs)
-        self.notify_observers()
-        return r
-
-    def __irshift__(self, *args, **kwargs):
-        r = np.ndarray.__irshift__(self, *args, **kwargs)
-        self.notify_observers()
-        return r
-
-
-    def __ixor__(self, *args, **kwargs):
-        r = np.ndarray.__ixor__(self, *args, **kwargs)
-        self.notify_observers()
-        return r
-
-
-    def __ipow__(self, *args, **kwargs):
-        r = np.ndarray.__ipow__(self, *args, **kwargs)
-        self.notify_observers()
-        return r
-
-
-    def __ifloordiv__(self, *args, **kwargs):
-        r = np.ndarray.__ifloordiv__(self, *args, **kwargs)
-        self.notify_observers()
-        return r
-
-
-    def __isub__(self, *args, **kwargs):
-        r = np.ndarray.__isub__(self, *args, **kwargs)
-        self.notify_observers()
-        return r
-
-
-    def __ior__(self, *args, **kwargs):
-        r = np.ndarray.__ior__(self, *args, **kwargs)
-        self.notify_observers()
-        return r
-
-
-    def __itruediv__(self, *args, **kwargs):
-        r = np.ndarray.__itruediv__(self, *args, **kwargs)
-        self.notify_observers()
-        return r
-
-
-    def __idiv__(self, *args, **kwargs):
-        r = np.ndarray.__idiv__(self, *args, **kwargs)
-        self.notify_observers()
-        return r
-
-
-    def __iand__(self, *args, **kwargs):
-        r = np.ndarray.__iand__(self, *args, **kwargs)
-        self.notify_observers()
-        return r
-
-
-    def __imod__(self, *args, **kwargs):
-        r = np.ndarray.__imod__(self, *args, **kwargs)
-        self.notify_observers()
-        return r
-
-
-    def __iadd__(self, *args, **kwargs):
-        r = np.ndarray.__iadd__(self, *args, **kwargs)
-        self.notify_observers()
-        return r
-
-
-    def __imul__(self, *args, **kwargs):
-        r = np.ndarray.__imul__(self, *args, **kwargs)
-        self.notify_observers()
-        return r
\ No newline at end of file
diff --git a/GPy/core/parameterization/param.py b/GPy/core/parameterization/param.py
index 60bdfe9d..182af902 100644
--- a/GPy/core/parameterization/param.py
+++ b/GPy/core/parameterization/param.py
@@ -4,7 +4,7 @@
 import itertools
 import numpy
 from parameter_core import OptimizationHandlable, adjust_name_for_printing
-from observable_array import ObsAr
+from array_core import ObsAr
 
 ###### printing
 __constraints_name__ = "Constraint"
diff --git a/GPy/kern/__init__.py b/GPy/kern/__init__.py
index a04c16a5..ef99e9a6 100644
--- a/GPy/kern/__init__.py
+++ b/GPy/kern/__init__.py
@@ -21,4 +21,5 @@ except ImportError:
 if sympy_available:
     from _src.symbolic import Symbolic
     from _src.heat_eqinit import Heat_eqinit
+    from _src.ode1_eq_lfm import Ode1_eq_lfm
     
diff --git a/GPy/kern/_src/symbolic.py b/GPy/kern/_src/symbolic.py
index 79abd91d..91193c5d 100644
--- a/GPy/kern/_src/symbolic.py
+++ b/GPy/kern/_src/symbolic.py
@@ -10,7 +10,7 @@ except ImportError:
 import numpy as np
 from kern import Kern
 from scipy.special import gammaln, gamma, erf, erfc, erfcx, polygamma
-from GPy.util.functions import normcdf, normcdfln, logistic, logisticln
+from GPy.util.functions import normcdf, normcdfln, logistic, logisticln, differfln
 from ...core.parameterization import Param
 
 class Symbolic(Kern):
@@ -39,6 +39,7 @@ class Symbolic(Kern):
                                'erf':erf, 'erfc':erfc,
                                'erfcx':erfcx,
                                'polygamma':polygamma,
+                               'differfln':differfln,
                                'normcdf':normcdf,
                                'normcdfln':normcdfln,
                                'logistic':logistic,
diff --git a/GPy/likelihoods/symbolic.py b/GPy/likelihoods/symbolic.py
index 65c90183..370f4599 100644
--- a/GPy/likelihoods/symbolic.py
+++ b/GPy/likelihoods/symbolic.py
@@ -117,21 +117,30 @@ if sympy_available:
             # generate the code for the pdf and derivatives
             self._gen_code()
 
+        def list_functions(self):
+            """Return a list of all symbolic functions in the model and their names."""
         def _gen_code(self):
             """Generate the code from the symbolic parts that will be used for likleihod computation."""
             # TODO: Check here whether theano is available and set up
             # functions accordingly.
-            self._log_pdf_function = lambdify(self.arg_list, self._sym_log_pdf, self.func_modules)
-
-            # compute code for derivatives 
-            self._derivative_code = {key: lambdify(self.arg_list, self._log_pdf_derivatives[key], self.func_modules) for key in self._log_pdf_derivatives.keys()}
-            self._second_derivative_code = {key: lambdify(self.arg_list, self._log_pdf_second_derivatives[key], self.func_modules) for key in self._log_pdf_second_derivatives.keys()}
-            self._third_derivative_code = {key: lambdify(self.arg_list, self._log_pdf_third_derivatives[key], self.func_modules) for key in self._log_pdf_third_derivatives.keys()}
-
+            symbolic_functions = [self._sym_log_pdf]
+            deriv_list = [self._log_pdf_derivatives, self._log_pdf_second_derivatives, self._log_pdf_third_derivatives]
+            symbolic_functions += [deriv[key] for key in sorted(deriv.keys()) for deriv in deriv_list]
             if self.missing_data:
-                self._missing_derivative_code = {key: lambdify(self.arg_list, self._missing_log_pdf_derivatives[key], self.func_modules) for key in self._missing_log_pdf_derivatives.keys()}
-                self._missing_second_derivative_code = {key: lambdify(self.arg_list, self._missing_log_pdf_second_derivatives[key], self.func_modules) for key in self._missing_log_pdf_second_derivatives.keys()}
-                self._missing_third_derivative_code = {key: lambdify(self.arg_list, self._missing_log_pdf_third_derivatives[key], self.func_modules) for key in self._missing_log_pdf_third_derivatives.keys()}
+                symbolic_functions+=[self._sym_missing_log_pdf]
+                deriv_list = [self._missing_log_pdf_derivatives, self._missing_log_pdf_second_derivatives, self._missing_log_pdf_third_derivatives]
+                symbolic_functions += [deriv[key] for key in sorted(deriv.keys()) for deriv in deriv_list]
+            # self._log_pdf_function = lambdify(self.arg_list, self._sym_log_pdf, self.func_modules)
+
+            # # compute code for derivatives 
+            # self._derivative_code = {key: lambdify(self.arg_list, self._log_pdf_derivatives[key], self.func_modules) for key in self._log_pdf_derivatives.keys()}
+            # self._second_derivative_code = {key: lambdify(self.arg_list, self._log_pdf_second_derivatives[key], self.func_modules) for key in self._log_pdf_second_derivatives.keys()}
+            # self._third_derivative_code = {key: lambdify(self.arg_list, self._log_pdf_third_derivatives[key], self.func_modules) for key in self._log_pdf_third_derivatives.keys()}
+
+            # if self.missing_data:
+            #     self._missing_derivative_code = {key: lambdify(self.arg_list, self._missing_log_pdf_derivatives[key], self.func_modules) for key in self._missing_log_pdf_derivatives.keys()}
+            #     self._missing_second_derivative_code = {key: lambdify(self.arg_list, self._missing_log_pdf_second_derivatives[key], self.func_modules) for key in self._missing_log_pdf_second_derivatives.keys()}
+            #     self._missing_third_derivative_code = {key: lambdify(self.arg_list, self._missing_log_pdf_third_derivatives[key], self.func_modules) for key in self._missing_log_pdf_third_derivatives.keys()}
 
             # TODO: compute EP code parts based on logZ. We need dlogZ/dmu, d2logZ/dmu2 and dlogZ/dtheta
 
diff --git a/GPy/util/functions.py b/GPy/util/functions.py
index 1cb1f7cb..3278182f 100644
--- a/GPy/util/functions.py
+++ b/GPy/util/functions.py
@@ -19,3 +19,7 @@ def normcdfln(x):
 
 def clip_exp(x):
     return np.where(x<lim_val, np.where(x>-lim_val, np.exp(x), epsilon), 1/epsilon)
+
+def differfln(x0, x1):
+    # this is a, hopefully!, a numerically more stable variant of log(erf(x0)-erf(x1)) = log(erfc(x1)-erfc(x0)).
+    return np.where(x0>x1, -x1*x1 + np.log(erfcx(x1)-np.exp(-x0**2+x1**2)*erfcx(x0)), -x0*x0 + np.log(np.exp(-x1**2+x0**2)*erfcx(x1) - erfcx(x0)))
diff --git a/GPy/util/symbolic.py b/GPy/util/symbolic.py
index 31aa1375..2cbff28f 100644
--- a/GPy/util/symbolic.py
+++ b/GPy/util/symbolic.py
@@ -1,10 +1,52 @@
 import sympy as sym
-from sympy import Function, S, oo, I, cos, sin, asin, log, erf, pi, exp, sqrt, sign, gamma,polygamma
-
+from sympy import Function, S, oo, I, cos, sin, asin, log, erf, pi, exp, sqrt, sign, gamma, polygamma
+from sympy.utilities.lambdify import lambdastr
+from sympy.utilities.iterables import numbered_symbols
 def stabilise(e):
     """Attempt to find the most numerically stable form of an expression"""
     return e #sym.expand(e)
 
+
+def gen_code(expressions, prefix='sub'):
+    """Generate code for the list of expressions provided using the common sub-expression eliminator."""
+    # First convert the expressions to a list.
+    # We should find the following type of expressions: 'function', 'derivative', 'second_derivative', 'third_derivative'. 
+
+    expression_list = []
+    expression_keys = []
+    function_code = {}
+    for key in expressions.key():
+        if key == 'function':
+            expression_list.append(expressions[key])
+            expression_keys.append(key)
+            function_code[key] = ''
+        elif key[-9:] == 'derivative':
+            function_code[key] = {}
+            for dkey in expressions[key]:
+                expression_list.append(expressions[key][dkey])
+                expression_keys.append([key, dkey])
+                function_code[key][dkey] = ''
+
+    symbols, functions = sym.cse(expression_list, numbered_symbols(prefix=prefix))
+    
+    # Create strings that lambda strings from the expressions. 
+    sub_expressions = []
+    for sub, expr in symbols:
+        arg_list = [e for e in expr.atoms() if e.is_Symbol]
+        sub_code += [lambdastr(sorted(arg_list), expr)]
+
+    function_expressions = []
+    for expr, keys in zip(functions, expression_keys):
+        arg_list = [e for e in expr.atoms() if e.is_Symbol]
+        function_code += [lambdastr(sorted(arg_list), expr)]
+
+    for key in function_code.key():
+        if key == 'function':
+            function_code[key] = 
+
+    return sub_code, func_code
+    
+
 class logistic(Function):
     """The logistic function as a symbolic function."""
     nargs = 1
@@ -123,23 +165,23 @@ class normal(Function):
     def _eval_is_real(self):
         return True
 
-class ln_diff_erf(Function):
+class differfln(Function):
     nargs = 2
 
     def fdiff(self, argindex=2):
         if argindex == 2:
             x0, x1 = self.args
-            return -2*exp(-x1**2)/(sqrt(pi)*(erf(x0)-erf(x1)))
+            return -2/(sqrt(pi)*(erfcx(x1)-exp(x1**2-x0**2)*erfcx(x0)))
         elif argindex == 1:
             x0, x1 = self.args
-            return 2.*exp(-x0**2)/(sqrt(pi)*(erf(x0)-erf(x1)))
+            return 2/(sqrt(pi)*(exp(x0**2-x1**2)*erfcx(x1)-erfcx(x0)))
         else:
             raise ArgumentIndexError(self, argindex)
         
     @classmethod
     def eval(cls, x0, x1):
         if x0.is_Number and x1.is_Number:            
-            return log(erf(x0)-erf(x1))
+            return log(erfc(x1)-erfc(x0))
 
 class dh_dd_i(Function):
     nargs = 5
@@ -304,7 +346,6 @@ class h(Function):
     
     @classmethod
     def eval(cls, t, tprime, d_i, d_j, l):
-        # putting in the is_Number stuff forces it to look for a fdiff method for derivative. If it's left out, then when asking for self.diff, it just does the diff on the eval symbolic terms directly. We want to avoid that because we are looking to ensure everything is numerically stable. Maybe it's because of the if statement that this happens? 
         if (t.is_Number
             and tprime.is_Number
             and d_i.is_Number