diff --git a/MathLib/ODE/Handles.h b/MathLib/ODE/Handles.h
new file mode 100644
index 0000000000000000000000000000000000000000..09920219630c21e9e2d61959481c661b2ebee273
--- /dev/null
+++ b/MathLib/ODE/Handles.h
@@ -0,0 +1,103 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef MATHLIB_ODE_HANDLES_H
+#define MATHLIB_ODE_HANDLES_H
+
+namespace MathLib
+{
+namespace detail
+{
+/// Function handles for N equations and arbitrary arguments.
+template <unsigned N, typename... FunctionArguments>
+struct Handles;
+
+/// Function handles for N equations and single argument.
+template <unsigned N, typename FunctionArgument>
+struct Handles<N, FunctionArgument> : public MathLib::FunctionHandles
+{
+	using Function = MathLib::Function<N, FunctionArgument>;
+	using JacobianFunction = MathLib::JacobianFunction<N, FunctionArgument>;
+
+	Handles(Function& f, JacobianFunction& df, FunctionArgument& arg)
+	    : f(f), df(df), _data(arg)
+	{
+	}
+
+	bool call(const double t, const double* const y,
+	          double* const ydot) override
+	{
+		// looks like f and df could be any callable object with suitable
+		// signature
+		// consider omission of data pointer and switch to std::function or
+		// alike
+		if (f)
+			return f(t, MappedConstVector<N>{y}, MappedVector<N>{ydot}, _data);
+		return true;
+	}
+
+	bool callJacobian(const double t, const double* const y, double* const ydot,
+	                  double* const jac) override
+	{
+		if (df)
+			return df(t,
+			          MappedConstVector<N>{y},
+			          MappedVector<N>{ydot},
+			          MappedMatrix<N, N>{jac},
+			          _data);
+		return true;
+	}
+
+	bool hasJacobian() const override { return df != nullptr; }
+	unsigned getNumEquations() const override { return N; }
+private:
+	Function f = nullptr;
+	JacobianFunction df = nullptr;
+	FunctionArgument& _data;
+};
+
+/// Function handles for N equations and no arguments.
+template <unsigned N>
+struct Handles<N> : public MathLib::FunctionHandles
+{
+	using Function = MathLib::Function<N>;
+	using JacobianFunction = MathLib::JacobianFunction<N>;
+
+	Handles(Function& f, JacobianFunction& df) : f(f), df(df) {}
+	bool call(const double t, const double* const y,
+	          double* const ydot) override
+	{
+		if (f)
+		{
+			return f(t, MappedConstVector<N>{y}, MappedVector<N>{ydot});
+		}
+		return true;
+	}
+
+	bool callJacobian(const double t, const double* const y, double* const ydot,
+	                  double* const jac) override
+	{
+		if (df)
+			return df(t,
+			          MappedConstVector<N>{y},
+			          MappedVector<N>{ydot},
+			          MappedMatrix<N, N>{jac});
+		return true;
+	}
+
+	bool hasJacobian() const override { return df != nullptr; }
+	unsigned getNumEquations() const override { return N; }
+	Function f = nullptr;
+	JacobianFunction df = nullptr;
+};
+
+}  // namespace detail
+}  // namespace MathLib
+
+#endif  // MATHLIB_ODE_HANDLES_H
diff --git a/MathLib/ODE/OdeSolverFactory.h b/MathLib/ODE/OdeSolverFactory.h
index 435d8326003e92e7f4818e61b66a1d8023545373..94ab4cd23ed45591c45303acd21fcc36577c478d 100644
--- a/MathLib/ODE/OdeSolverFactory.h
+++ b/MathLib/ODE/OdeSolverFactory.h
@@ -15,97 +15,12 @@
 #include "BaseLib/ConfigTree.h"
 
 #include "OdeSolver.h"
+#include "Handles.h"
 
 #include "CVodeSolver.h"
 
 namespace MathLib
 {
-namespace detail
-{
-template <unsigned N, typename... FunctionArguments>
-struct Handles;
-
-template <unsigned N, typename FunctionArgument>
-struct Handles<N, FunctionArgument> : public MathLib::FunctionHandles
-{
-	using Function = MathLib::Function<N, FunctionArgument>;
-	using JacobianFunction = MathLib::JacobianFunction<N, FunctionArgument>;
-
-	Handles(Function& f, JacobianFunction& df, FunctionArgument& arg)
-	    : f(f), df(df), _data(arg)
-	{
-	}
-
-	bool call(const double t, const double* const y,
-	          double* const ydot) override
-	{
-		// looks like f and df could be any callable object with suitable
-		// signature
-		// consider omission of data pointer and switch to std::function or
-		// alike
-		if (f)
-			return f(t, MappedConstVector<N>{y}, MappedVector<N>{ydot}, _data);
-		return true;
-	}
-
-	bool callJacobian(const double t, const double* const y, double* const ydot,
-	                  double* const jac) override
-	{
-		if (df)
-			return df(t,
-			          MappedConstVector<N>{y},
-			          MappedVector<N>{ydot},
-			          MappedMatrix<N, N>{jac /*, order*/},
-			          _data);
-		return true;
-	}
-
-	bool hasJacobian() const override { return df != nullptr; }
-	unsigned getNumEquations() const override { return N; }
-private:
-	Function f = nullptr;
-	JacobianFunction df = nullptr;
-	FunctionArgument& _data;
-};
-
-template <unsigned N>
-struct Handles<N> : public MathLib::FunctionHandles
-{
-	using Function = MathLib::Function<N>;
-	using JacobianFunction = MathLib::JacobianFunction<N>;
-
-	Handles(Function& f, JacobianFunction& df) : f(f), df(df) {}
-	bool call(const double t, const double* const y,
-	          double* const ydot) override
-	{
-		if (f)
-		{
-			// auto ydot_ = Eigen::Map<Eigen::Matrix<double, N, 1>>{y};
-			// auto ydot_ = Eigen::Map<Eigen::Matrix<double, N, 1>>{ydot};
-			return f(t, MappedConstVector<N>{y}, MappedVector<N>{ydot});
-		}
-		return true;
-	}
-
-	bool callJacobian(const double t, const double* const y, double* const ydot,
-	                  double* const jac) override
-	{
-		if (df)
-			return df(t,
-			          MappedConstVector<N>{y},
-			          MappedVector<N>{ydot},
-			          MappedMatrix<N, N>{jac /*, order*/});
-		return true;
-	}
-
-	bool hasJacobian() const override { return df != nullptr; }
-	unsigned getNumEquations() const override { return N; }
-	Function f = nullptr;
-	JacobianFunction df = nullptr;
-};
-
-}  // namespace detail
-
 template <unsigned NumEquations, typename... FunctionArguments>
 std::unique_ptr<OdeSolver<NumEquations, FunctionArguments...>> createOdeSolver(
     BaseLib::ConfigTree const& config);