diff --git a/MathLib/CMakeLists.txt b/MathLib/CMakeLists.txt
index a97ed510448c9db74449c9ee6a996f16f6ff5daf..8a8b224fcea564e75d9bfa73b76a702ed20b6a17 100644
--- a/MathLib/CMakeLists.txt
+++ b/MathLib/CMakeLists.txt
@@ -23,6 +23,9 @@ set(SOURCES ${SOURCES} ${SOURCES_LINALG_SOLVERS})
 GET_SOURCE_FILES(SOURCES_LINALG_PRECOND LinAlg/Preconditioner)
 set(SOURCES ${SOURCES} ${SOURCES_LINALG_PRECOND})
 
+GET_SOURCE_FILES(SOURCES_ODE ODE)
+set(SOURCES ${SOURCES} ${SOURCES_ODE})
+
 if(OGS_USE_EIGEN)
     GET_SOURCE_FILES(SOURCES_LINALG_EIGEN LinAlg/Eigen)
     set(SOURCES ${SOURCES} ${SOURCES_LINALG_EIGEN})
@@ -62,6 +65,8 @@ set_target_properties(MathLib PROPERTIES LINKER_LANGUAGE CXX)
 
 target_link_libraries(MathLib
     AssemblerLib
+    sundials_cvode
+    sundials_nvecserial
 )
 
 if(METIS_FOUND)
diff --git a/MathLib/ODE/CVodeSolver.cpp b/MathLib/ODE/CVodeSolver.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..aef35994d15d5df915d122691b3d6da633b970e5
--- /dev/null
+++ b/MathLib/ODE/CVodeSolver.cpp
@@ -0,0 +1,363 @@
+#include "CVodeSolver.h"
+
+#include <cassert>
+#include <logog/include/logog.hpp>
+
+#include <cvode/cvode.h>             /* prototypes for CVODE fcts., consts. */
+#include <nvector/nvector_serial.h>  /* serial N_Vector types, fcts., macros */
+#include <cvode/cvode_dense.h>       /* prototype for CVDense */
+#include <sundials/sundials_dense.h> /* definitions DlsMat DENSE_ELEM */
+#include <sundials/sundials_types.h> /* definition of type realtype */
+
+
+namespace
+{
+void printStats(void* cvode_mem)
+{
+	long int nst, nfe, nsetups, nje, nfeLS, nni, ncfn, netf, nge;
+	int flag;
+
+	flag = CVodeGetNumSteps(cvode_mem, &nst);
+	// check_flag(&flag, "CVodeGetNumSteps", 1);
+	flag = CVodeGetNumRhsEvals(cvode_mem, &nfe);
+	// check_flag(&flag, "CVodeGetNumRhsEvals", 1);
+	flag = CVodeGetNumLinSolvSetups(cvode_mem, &nsetups);
+	// check_flag(&flag, "CVodeGetNumLinSolvSetups", 1);
+	flag = CVodeGetNumErrTestFails(cvode_mem, &netf);
+	// check_flag(&flag, "CVodeGetNumErrTestFails", 1);
+	flag = CVodeGetNumNonlinSolvIters(cvode_mem, &nni);
+	// check_flag(&flag, "CVodeGetNumNonlinSolvIters", 1);
+	flag = CVodeGetNumNonlinSolvConvFails(cvode_mem, &ncfn);
+	// check_flag(&flag, "CVodeGetNumNonlinSolvConvFails", 1);
+
+	flag = CVDlsGetNumJacEvals(cvode_mem, &nje);
+	// check_flag(&flag, "CVDlsGetNumJacEvals", 1);
+	flag = CVDlsGetNumRhsEvals(cvode_mem, &nfeLS);
+	// check_flag(&flag, "CVDlsGetNumRhsEvals", 1);
+
+	flag = CVodeGetNumGEvals(cvode_mem, &nge);
+	// check_flag(&flag, "CVodeGetNumGEvals", 1);
+
+	DBUG("\nFinal Statistics:");
+	DBUG("nst = %-6ld nfe  = %-6ld nsetups = %-6ld nfeLS = %-6ld nje = %ld",
+	     nst, nfe, nsetups, nfeLS, nje);
+	DBUG("nni = %-6ld ncfn = %-6ld netf = %-6ld nge = %ld\n", nni, ncfn, netf,
+	     nge);
+}
+}
+
+namespace MathLib
+{
+/**
+ * This class provides concrete access to Sundials' CVode solver.
+ *
+ * OdeSolver (implicitly bounds checked, agnostic to concrete ODE solver)
+ *  |
+ *  | Dynamic polymorphism
+ *  v
+ * ConcreteOdeSolver (implicitly bounds checked, interfaces with a specific
+ * library)
+ *  |
+ *  | Forward calls, disable bounds checking, no need for templates anymore
+ *  v
+ * Implementation = CVodeSolverInternal (no templates)
+ *  |
+ *  | Pimpl (hide implementation, do not include 3rd party libs in header)
+ *  v
+ * CVodeSolverImpl
+ *
+ * This scheme might be a general way for accessing 3rd party libraries.
+ */
+class CVodeSolverImpl
+{
+	static_assert(std::is_same<realtype, double>::value,
+	              "cvode's realtype is not the same as double");
+
+public:
+	CVodeSolverImpl(BaseLib::ConfigTree const& config);
+	void init(const unsigned num_equations);
+
+	friend class CVodeSolver;
+	~CVodeSolverImpl();
+
+private:
+	void setFunction(FunctionHandles* f);
+
+	void preSolve();
+	void solve(const double t_end);
+
+	double const* getSolution() const { return NV_DATA_S(_y); }
+	double getTime() const { return _t; }
+	bool getYDot(const double t, double const* const y, double* const ydot);
+	void setTolerance(const double* abstol, const double reltol);
+	void setTolerance(const double abstol, const double reltol);
+	void setIC(const double t0, double const* const y0);
+
+private:
+	N_Vector _y = nullptr;
+
+	realtype _t;
+
+	N_Vector _abstol = nullptr;
+	realtype _reltol;
+
+	unsigned _num_equations;
+	void* _cvode_mem;
+
+	FunctionHandles* _f;
+
+	int _linear_multistep_method = CV_ADAMS;
+	int _nonlinear_solver_iteration = CV_FUNCTIONAL;
+};
+
+CVodeSolverImpl::CVodeSolverImpl(const BaseLib::ConfigTree& config)
+{
+	if (auto const param =
+	        config.getConfParamOptional<std::string>("linear_multistep_method"))
+	{
+		DBUG("setting linear multistep method (config: %s)", param->c_str());
+
+		if (*param == "Adams")
+		{
+			_linear_multistep_method = CV_ADAMS;
+		}
+		else if (*param == "BDF")
+		{
+			_linear_multistep_method = CV_BDF;
+		}
+		else
+		{
+			ERR("unknown linear multistep method: %s", param->c_str());
+			std::abort();
+		}
+	}
+
+	if (auto const param = config.getConfParamOptional<std::string>(
+	        "nonlinear_solver_iteration"))
+	{
+		DBUG("setting nonlinear solver iteration (config: %s)", param->c_str());
+
+		if (*param == "Functional")
+		{
+			_nonlinear_solver_iteration = CV_FUNCTIONAL;
+		}
+		else if (*param == "Newton")
+		{
+			_nonlinear_solver_iteration = CV_NEWTON;
+		}
+		else
+		{
+			ERR("unknown nonlinear solver iteration: %s", param->c_str());
+			std::abort();
+		}
+	}
+}
+
+void CVodeSolverImpl::init(const unsigned num_equations)
+{
+	_y = N_VNew_Serial(num_equations);
+	_abstol = N_VNew_Serial(num_equations);
+	_num_equations = num_equations;
+
+	_cvode_mem =
+	    CVodeCreate(_linear_multistep_method, _nonlinear_solver_iteration);
+
+	assert(_cvode_mem != nullptr && _y != nullptr && _abstol != nullptr);
+}
+
+void CVodeSolverImpl::setTolerance(const double* abstol, const double reltol)
+{
+	for (unsigned i = 0; i < _num_equations; ++i)
+	{
+		NV_Ith_S(_abstol, i) = abstol[i];
+	}
+
+	_reltol = reltol;
+}
+
+void CVodeSolverImpl::setTolerance(const double abstol, const double reltol)
+{
+	for (unsigned i = 0; i < _num_equations; ++i)
+	{
+		NV_Ith_S(_abstol, i) = abstol;
+	}
+
+	_reltol = reltol;
+}
+
+void CVodeSolverImpl::setFunction(FunctionHandles* f)
+{
+	_f = f;
+	assert(_num_equations == f->getNumEquations());
+
+	auto f_wrapped = [](const realtype t, const N_Vector y, N_Vector ydot,
+	                    void* function_handles) -> int
+	{
+		bool successful = static_cast<FunctionHandles*>(function_handles)
+		                      ->call(t, NV_DATA_S(y), NV_DATA_S(ydot));
+		return successful ? 0 : 1;
+	};
+
+	// TODO: check not run twice! move this call somewhere else
+	int flag = CVodeInit(_cvode_mem, f_wrapped, 0.0, _y);
+	(void)flag;
+}
+
+void CVodeSolverImpl::setIC(const double t0, double const* const y0)
+{
+	for (unsigned i = 0; i < _num_equations; ++i)
+	{
+		NV_Ith_S(_y, i) = y0[i];
+	}
+
+	_t = t0;
+}
+
+void CVodeSolverImpl::preSolve()
+{
+	assert(_f != nullptr && "ode function handle was not provided");
+
+	// int flag = CVodeInit(_cvode_mem, f_wrapped, _t, _y); // TODO: consider
+	// CVodeReInit()!
+	int flag =
+	    CVodeReInit(_cvode_mem, _t, _y);  // TODO: consider CVodeReInit()!
+	// if (check_flag(&flag, "CVodeInit", 1)) return(1);
+
+	flag = CVodeSetUserData(_cvode_mem, static_cast<void*>(_f));
+
+	/* Call CVodeSVtolerances to specify the scalar relative tolerance
+	 * and vector absolute tolerances */
+	flag = CVodeSVtolerances(_cvode_mem, _reltol, _abstol);
+	// if (check_flag(&flag, "CVodeSVtolerances", 1)) return(1);
+
+	/* Call CVDense to specify the CVDENSE dense linear solver */
+	flag = CVDense(_cvode_mem, _num_equations);
+	// if (check_flag(&flag, "CVDense", 1)) return(1);
+
+	if (_f->hasJacobian())
+	{
+		auto df_wrapped = [](
+		    const long /*N*/, const realtype t, const N_Vector y,
+		    const N_Vector ydot, const DlsMat jac, void* function_handles,
+		    N_Vector /*tmp1*/, N_Vector /*tmp2*/, N_Vector /*tmp3*/
+		    ) -> int
+		{
+			// Caution: by calling the DENSE_COL() macro we assume that matrices
+			//          are stored contiguously in memory!
+			bool successful =
+			    static_cast<FunctionHandles*>(function_handles)
+			        ->callJacobian(t, NV_DATA_S(y), NV_DATA_S(ydot),
+			                       DENSE_COL(jac, 0),
+			                       BaseLib::StorageOrder::ColumnMajor);
+			return successful ? 0 : 1;
+		};
+
+		flag = CVDlsSetDenseJacFn(_cvode_mem, df_wrapped);
+		// if (check_flag(&flag, "CVDlsSetDenseJacFn", 1)) return(1);
+	}
+}
+
+void CVodeSolverImpl::solve(const double t_end)
+{
+	realtype t_reached;
+	int flag = CVode(_cvode_mem, t_end, _y, &t_reached, CV_NORMAL);
+	_t = t_reached;
+	// std::cout << "result at time " << t << " is " << NV_Ith_S(y,0) <<
+	// std::endl;
+	if (flag != CV_SUCCESS)
+	{
+		// std::cerr << "ERROR at " << __FUNCTION__ << ":" << __LINE__ <<
+		// std::endl;
+	}
+}
+
+bool CVodeSolverImpl::getYDot(const double t, double const* const y,
+                              double* const ydot)
+{
+	if (_f != nullptr)
+	{
+		return _f->call(t, y, ydot);
+	}
+
+	return false;
+}
+
+CVodeSolverImpl::~CVodeSolverImpl()
+{
+	printStats(_cvode_mem);
+
+	if (_y)
+	{
+		N_VDestroy_Serial(_y);
+		N_VDestroy_Serial(_abstol);
+	}
+
+	if (_cvode_mem)
+	{
+		CVodeFree(&_cvode_mem);
+	}
+}
+
+CVodeSolverInternal::CVodeSolverInternal(BaseLib::ConfigTree const& config)
+    : _impl(new CVodeSolverImpl(config))
+{
+}
+
+void CVodeSolverInternal::init(const unsigned num_equations)
+{
+	_impl->init(num_equations);
+}
+
+void CVodeSolverInternal::setTolerance(const double* abstol,
+                                       const double reltol)
+{
+	_impl->setTolerance(abstol, reltol);
+}
+
+void CVodeSolverInternal::setTolerance(const double abstol, const double reltol)
+{
+	_impl->setTolerance(abstol, reltol);
+}
+
+void CVodeSolverInternal::setFunction(FunctionHandles* f)
+{
+	_impl->setFunction(f);
+}
+
+void CVodeSolverInternal::setIC(const double t0, double const* const y0)
+{
+	_impl->setIC(t0, y0);
+}
+
+void CVodeSolverInternal::preSolve()
+{
+	_impl->preSolve();
+}
+
+void CVodeSolverInternal::solve(const double t_end)
+{
+	_impl->solve(t_end);
+}
+
+double const* CVodeSolverInternal::getSolution() const
+{
+	return _impl->getSolution();
+}
+
+bool CVodeSolverInternal::getYDot(const double t, double const* const y,
+                                  double* const ydot) const
+{
+	return _impl->getYDot(t, y, ydot);
+}
+
+double CVodeSolverInternal::getTime() const
+{
+	return _impl->getTime();
+}
+
+CVodeSolverInternal::~CVodeSolverInternal()
+{
+	delete _impl;
+}
+
+}  // namespace MathLib
diff --git a/MathLib/ODE/CVodeSolver.h b/MathLib/ODE/CVodeSolver.h
new file mode 100644
index 0000000000000000000000000000000000000000..a1309ee82ab0c07db50c6c583702e2dc0fad25aa
--- /dev/null
+++ b/MathLib/ODE/CVodeSolver.h
@@ -0,0 +1,47 @@
+#ifndef MATHLIB_CVODESOLVER_H
+#define MATHLIB_CVODESOLVER_H
+
+#include "BaseLib/ConfigTree.h"
+
+#include "declarations.h"
+
+namespace MathLib
+{
+class CVodeSolverImpl;
+
+/**
+ * ODE solver, general, pointer based implementation. No implicit bounds
+ * checking
+ *
+ * For internal use only.
+ */
+class CVodeSolverInternal
+{
+protected:
+	CVodeSolverInternal(BaseLib::ConfigTree const& config);
+	void init(const unsigned num_equations);
+
+	void setTolerance(double const* const abstol, const double reltol);
+	void setTolerance(const double abstol, const double reltol);
+
+	void setFunction(FunctionHandles* f);
+
+	void setIC(const double t0, double const* const y0);
+
+	void preSolve();
+	void solve(const double t_end);
+
+	double const* getSolution() const;
+	double getTime() const;
+	bool getYDot(const double t, double const* const y,
+	             double* const ydot) const;
+
+	~CVodeSolverInternal();
+
+private:
+	CVodeSolverImpl* _impl;  ///< pimpl idiom hides implementation
+};
+
+}  // namespace MathLib
+
+#endif  // MATHLIB_CVODESOLVER_H
diff --git a/MathLib/ODE/OdeSolver.h b/MathLib/ODE/OdeSolver.h
new file mode 100644
index 0000000000000000000000000000000000000000..ad79a1bd3014467d7da0aafed5b8cd88138d0a9e
--- /dev/null
+++ b/MathLib/ODE/OdeSolver.h
@@ -0,0 +1,50 @@
+#ifndef MATHLIB_ODESOLVER_H
+#define MATHLIB_ODESOLVER_H
+
+#include <array>
+
+#include "declarations.h"
+
+namespace MathLib
+{
+/**
+ * ODE solver Interface.
+ *
+ * This class provides an abstract interface for ODE solvers.
+ * It provides type-safe and array-bounds checked access to external
+ * ODE solver libraries. However, it is agnostic to the specific solver used.
+ */
+template <unsigned NumEquations, typename... FunctionArguments>
+class OdeSolver
+{
+public:
+	using Arr = std::array<double, NumEquations>;
+	using ConstArrRef = BaseLib::ArrayRef<const double, NumEquations>;
+	using Function = MathLib::Function<NumEquations, FunctionArguments...>;
+	using JacobianFunction =
+	    MathLib::JacobianFunction<NumEquations, FunctionArguments...>;
+
+	virtual void init() = 0;
+
+	virtual void setTolerance(const Arr& abstol, const double reltol) = 0;
+	virtual void setTolerance(const double abstol, const double reltol) = 0;
+
+	virtual void setFunction(Function f, JacobianFunction df,
+	                         FunctionArguments*... args) = 0;
+
+	virtual void setIC(const double t0, const Arr& y0) = 0;
+
+	virtual void preSolve() = 0;
+	virtual void solve(const double t) = 0;
+
+	virtual unsigned getNumEquations() const { return NumEquations; }
+	virtual ConstArrRef getSolution() const = 0;
+	virtual double getTime() const = 0;
+	virtual Arr getYDot(const double t, const Arr& y) const = 0;
+
+	virtual ~OdeSolver() = default;
+};
+
+}  // namespace MathLib
+
+#endif  // MATHLIB_ODESOLVER_H
diff --git a/MathLib/ODE/OdeSolverFactory.h b/MathLib/ODE/OdeSolverFactory.h
new file mode 100644
index 0000000000000000000000000000000000000000..72964b296e29cbdb384261178b8eb10c66a22633
--- /dev/null
+++ b/MathLib/ODE/OdeSolverFactory.h
@@ -0,0 +1,204 @@
+#pragma once
+
+#include <memory>
+
+#include "BaseLib/ConfigTree.h"
+
+#include "OdeSolver.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>;
+
+	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,
+			         BaseLib::ArrayRef<const double, N>{y},
+			         BaseLib::ArrayRef<double, N>{ydot},
+			         *_data);
+		return true;
+	}
+
+	bool callJacobian(const double t, const double* const y,
+	                  const double* const ydot, double* const jac,
+	                  BaseLib::StorageOrder order) override
+	{
+		if (df)
+			return df(t,
+			          BaseLib::ArrayRef<const double, N>{y},
+			          BaseLib::ArrayRef<const double, N>{ydot},
+			          BaseLib::MatrixRef<double, N, N>{jac, order},
+			          *_data);
+		return true;
+	}
+
+	bool hasJacobian() const override { return df != nullptr; }
+	unsigned getNumEquations() const override { return N; }
+	void setArguments(FunctionArgument* arg)
+	{
+		assert(arg != nullptr);
+		_data = arg;
+	}
+
+	// TODO: make private
+	Function f = nullptr;
+	JacobianFunction df = nullptr;
+
+private:
+	FunctionArgument* _data = nullptr;
+};
+
+template <unsigned N>
+struct Handles<N> : public MathLib::FunctionHandles
+{
+	using Function = MathLib::Function<N>;
+	using JacobianFunction = MathLib::JacobianFunction<N>;
+
+	bool call(const double t, const double* const y,
+	          double* const ydot) override
+	{
+		if (f)
+			return f(t,
+			         BaseLib::ArrayRef<const double, N>{y},
+			         BaseLib::ArrayRef<double, N>{ydot});
+		return true;
+	}
+
+	bool callJacobian(const double t, const double* const y,
+	                  const double* const ydot, double* const jac,
+	                  BaseLib::StorageOrder order) override
+	{
+		if (df)
+			return df(t,
+			          BaseLib::ArrayRef<const double, N>{y},
+			          BaseLib::ArrayRef<const double, N>{ydot},
+			          BaseLib::MatrixRef<double, N, N>{jac, order});
+		return true;
+	}
+
+	bool hasJacobian() const override { return df != nullptr; }
+	unsigned getNumEquations() const override { return N; }
+	void setArguments() const {}
+	Function f = nullptr;
+	JacobianFunction df = nullptr;
+};
+
+}  // namespace detail
+
+template <unsigned NumEquations, typename... FunctionArguments>
+std::unique_ptr<OdeSolver<NumEquations, FunctionArguments...>> createOdeSolver(
+    BaseLib::ConfigTree const& config);
+
+/**
+ * ODE solver with a bounds-safe interface.
+ *
+ * This class makes contact between the abstract \c OdeSolver interface and a
+ * certain solver \c Implementation.
+ *
+ * The interface of this class inherits the array bounds checking from \c
+ * OdeSolver.
+ * Its methods forward calls to the \c Implementation erasing array bounds info
+ * by
+ * passing \c std::array as raw pointer.
+ *
+ * This way the \c Implementation does not need to be templated.
+ */
+template <unsigned NumEquations, typename Implementation,
+          typename... FunctionArguments>
+class ConcreteOdeSolver final
+    : public OdeSolver<NumEquations, FunctionArguments...>,
+      private Implementation
+{
+public:
+	using Interface = OdeSolver<NumEquations, FunctionArguments...>;
+	using Arr = typename Interface::Arr;
+	using ConstArrRef = typename Interface::ConstArrRef;
+	using Function = typename Interface::Function;
+	using JacobianFunction = typename Interface::JacobianFunction;
+
+	void init() override
+	{
+		Implementation::init(NumEquations);
+		Implementation::setFunction(&_handles);
+	}
+
+	void setTolerance(const Arr& abstol, const double reltol) override
+	{
+		Implementation::setTolerance(abstol.data(), reltol);
+	}
+
+	void setTolerance(const double abstol, const double reltol) override
+	{
+		Implementation::setTolerance(abstol, reltol);
+	}
+
+	void setFunction(Function f, JacobianFunction df,
+	                 FunctionArguments*... args) override
+	{
+		_handles.f = f;
+		_handles.df = df;
+		_handles.setArguments(args...);
+	}
+
+	void setIC(const double t0, const Arr& y0) override
+	{
+		Implementation::setIC(t0, y0.data());
+	}
+
+	void preSolve() { Implementation::preSolve(); }
+	void solve(const double t) override { Implementation::solve(t); }
+	ConstArrRef getSolution() const override
+	{
+		return ConstArrRef(Implementation::getSolution());
+	}
+
+	double getTime() const override { return Implementation::getTime(); }
+	Arr getYDot(const double t, const Arr& y) const override
+	{
+		Arr ydot;
+		Implementation::getYDot(t, y.data(), ydot.data());
+		return ydot;
+	}
+
+private:
+	/// instances of this class shall only be constructed by
+	/// the friend function listed below
+	ConcreteOdeSolver(BaseLib::ConfigTree const& config)
+	    : Implementation{config}
+	{
+	}
+
+	detail::Handles<NumEquations, FunctionArguments...> _handles;
+
+	friend std::unique_ptr<OdeSolver<NumEquations, FunctionArguments...>>
+	createOdeSolver<NumEquations, FunctionArguments...>(
+	    BaseLib::ConfigTree const& config);
+};
+
+template <unsigned NumEquations, typename... FunctionArguments>
+std::unique_ptr<OdeSolver<NumEquations, FunctionArguments...>> createOdeSolver(
+    BaseLib::ConfigTree const& config)
+{
+	return std::unique_ptr<OdeSolver<NumEquations, FunctionArguments...>>(
+	    new ConcreteOdeSolver<NumEquations, CVodeSolverInternal,
+	                          FunctionArguments...>(config));
+}
+
+}  // namespace MathLib
diff --git a/MathLib/ODE/declarations.h b/MathLib/ODE/declarations.h
new file mode 100644
index 0000000000000000000000000000000000000000..4799a37d0f64c5860a5e9f03f35149fa084c3804
--- /dev/null
+++ b/MathLib/ODE/declarations.h
@@ -0,0 +1,41 @@
+#pragma once
+
+#include <cassert>
+
+#include "BaseLib/ArrayRef.h"
+#include "BaseLib/MatrixRef.h"
+
+namespace MathLib
+{
+template <unsigned N, typename... FunctionArguments>
+using Function = bool (*)(const double t,
+                          BaseLib::ArrayRef<const double, N> y,
+                          BaseLib::ArrayRef<double, N> ydot,
+                          FunctionArguments&... arg);
+
+template <unsigned N, typename... FunctionArguments>
+using JacobianFunction = bool (*)(const double t,
+                                  BaseLib::ArrayRef<const double, N> y,
+                                  BaseLib::ArrayRef<const double, N> ydot,
+                                  BaseLib::MatrixRef<double, N, N> jac,
+                                  FunctionArguments&... arg);
+
+// This is an internal detail
+class FunctionHandles
+{
+public:
+	virtual bool call(const double t, double const* const y,
+	                  double* const ydot) = 0;
+	virtual bool callJacobian(const double t,
+	                          double const* const y,
+	                          double const* const ydot,
+	                          double* const jac,
+	                          BaseLib::StorageOrder order) = 0;
+
+	virtual bool hasJacobian() const = 0;
+
+	virtual unsigned getNumEquations() const = 0;
+
+	virtual ~FunctionHandles() = default;
+};
+}
diff --git a/Tests/MathLib/TestCVode.cpp b/Tests/MathLib/TestCVode.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..f691b13f6f8ba242af3013f939a5b514d22d1e84
--- /dev/null
+++ b/Tests/MathLib/TestCVode.cpp
@@ -0,0 +1,211 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/LICENSE.txt
+ */
+
+#include <gtest/gtest.h>
+
+#include "MathLib/ODE/CVodeSolver.h"
+#include "MathLib/ODE/OdeSolver.h"
+#include "MathLib/ODE/OdeSolverFactory.h"
+
+#include <cstdio>
+
+bool f(const double,
+       BaseLib::ArrayRef<const double, 1> y,
+       BaseLib::ArrayRef<double, 1> ydot)
+{
+	if (y[0] <= 0.0) return false;
+
+	ydot[0] = -15.0 * y[0];
+	return true;
+}
+
+bool df(const double /*t*/,
+        BaseLib::ArrayRef<const double, 1> y,
+        BaseLib::ArrayRef<const double, 1> /*ydot*/,
+        BaseLib::MatrixRef<double, 1, 1> jac)
+{
+	if (y[0] <= 0.0) return false;
+
+	jac(0, 0) = -15.0;
+	return true;
+}
+
+struct ExtraData
+{
+	double value = 15.0;
+};
+
+bool f_extra(const double,
+             BaseLib::ArrayRef<const double, 1> y,
+             BaseLib::ArrayRef<double, 1> ydot,
+             ExtraData& data)
+{
+	if (y[0] <= 0.0) return false;
+
+	ydot[0] = -data.value * y[0];
+	return true;
+}
+
+TEST(MathLibCVodeTest, Exponential)
+{
+	// initial values
+	const double y0 = 1.0;
+	const double t0 = 0.0;
+
+	BaseLib::ConfigTree config(boost::property_tree::ptree{}, "",
+	                           BaseLib::ConfigTree::onerror,
+	                           BaseLib::ConfigTree::onwarning);
+	auto ode_solver = MathLib::createOdeSolver<1>(config);
+
+	ode_solver->init();
+	ode_solver->setTolerance(1e-8, 1e-6);
+
+	ode_solver->setFunction(f, nullptr);
+
+	ode_solver->setIC(t0, {y0});
+
+	ode_solver->preSolve();
+
+	const double dt = 1e-1;
+
+	for (unsigned i = 1; i <= 10; ++i)
+	{
+		const double time = dt * i;
+
+		ode_solver->solve(time);
+
+		auto const y = ode_solver->getSolution();
+		double time_reached = ode_solver->getTime();
+
+		std::printf("t: %14.7g, y: %14.7g, diff: %14.7g\n", time_reached, y[0],
+		            y[0] - exp(-15.0 * time_reached));
+
+		ASSERT_NEAR(time, time_reached, std::numeric_limits<double>::epsilon());
+		// std::printf("time: %g\n", time_reached);
+	}
+}
+
+TEST(MathLibCVodeTest, ExponentialExtraData)
+{
+	// initial values
+	const double y0 = 1.0;
+	const double t0 = 0.0;
+
+	BaseLib::ConfigTree config(boost::property_tree::ptree{}, "",
+	                           BaseLib::ConfigTree::onerror,
+	                           BaseLib::ConfigTree::onwarning);
+	auto ode_solver = MathLib::createOdeSolver<1, ExtraData>(config);
+
+	ode_solver->init();
+	ode_solver->setTolerance(1e-8, 1e-6);
+
+	ExtraData data;
+	ode_solver->setFunction(f_extra, nullptr, &data);
+
+	ode_solver->setIC(t0, {y0});
+	ode_solver->preSolve();
+
+	const double dt = 1e-1;
+
+	for (unsigned i = 1; i <= 10; ++i)
+	{
+		const double time = dt * i;
+
+		ode_solver->solve(time);
+
+		auto const y = ode_solver->getSolution();
+		double time_reached = ode_solver->getTime();
+
+		std::printf("t: %14.7g, y: %14.7g, diff: %14.7g\n", time_reached, y[0],
+		            y[0] - exp(-15.0 * time_reached));
+
+		ASSERT_NEAR(time, time_reached, std::numeric_limits<double>::epsilon());
+		// std::printf("time: %g\n", time_reached);
+	}
+}
+
+TEST(MathLibCVodeTest, ExponentialWithJacobian)
+{
+	// initial values
+	const double y0 = 1.0;
+	const double t0 = 0.0;
+
+	BaseLib::ConfigTree config(boost::property_tree::ptree{}, "",
+	                           BaseLib::ConfigTree::onerror,
+	                           BaseLib::ConfigTree::onwarning);
+	auto ode_solver = MathLib::createOdeSolver<1>(config);
+
+	ode_solver->init();
+	ode_solver->setTolerance(1e-10, 1e-8);
+
+	ode_solver->setFunction(f, df);
+
+	ode_solver->setIC(t0, {y0});
+
+	ode_solver->preSolve();
+
+	const double dt = 1e-1;
+
+	for (unsigned i = 1; i <= 10; ++i)
+	{
+		const double time = dt * i;
+
+		ode_solver->solve(time);
+
+		auto const y = ode_solver->getSolution();
+		double time_reached = ode_solver->getTime();
+
+		std::printf("t: %14.7g, y: %14.7g, diff: %14.7g\n", time_reached, y[0],
+		            y[0] - exp(-15.0 * time_reached));
+
+		ASSERT_NEAR(time, time_reached, std::numeric_limits<double>::epsilon());
+		// std::printf("time: %g\n", time_reached);
+	}
+}
+
+TEST(MathLibCVodeTest, ExponentialWithJacobianNewton)
+{
+	// initial values
+	const double y0 = 1.0;
+	const double t0 = 0.0;
+
+	boost::property_tree::ptree conf;
+	conf.put("linear_multistep_method", "BDF");
+	conf.put("nonlinear_solver_iteration", "Newton");
+	BaseLib::ConfigTree config(conf, "", BaseLib::ConfigTree::onerror,
+	                           BaseLib::ConfigTree::onwarning);
+
+	auto ode_solver = MathLib::createOdeSolver<1>(config);
+
+	ode_solver->init();
+	ode_solver->setTolerance(1e-6, 1e-6);
+
+	ode_solver->setFunction(f, df);
+
+	ode_solver->setIC(t0, {y0});
+
+	ode_solver->preSolve();
+
+	const double dt = 1e-1;
+
+	for (unsigned i = 1; i <= 10; ++i)
+	{
+		const double time = dt * i;
+
+		ode_solver->solve(time);
+
+		auto const y = ode_solver->getSolution();
+		double time_reached = ode_solver->getTime();
+
+		std::printf("t: %14.7g, y: %14.7g, diff: %14.7g\n", time_reached, y[0],
+		            y[0] - exp(-15.0 * time_reached));
+
+		ASSERT_NEAR(time, time_reached, std::numeric_limits<double>::epsilon());
+		// std::printf("time: %g\n", time_reached);
+	}
+}