diff --git a/NumLib/ODESolver/NonlinearSolver-impl.h b/NumLib/ODESolver/NonlinearSolver-impl.h
new file mode 100644
index 0000000000000000000000000000000000000000..5fc9b8ed11add845be222d6f480e3a999a8bb6ef
--- /dev/null
+++ b/NumLib/ODESolver/NonlinearSolver-impl.h
@@ -0,0 +1,135 @@
+/**
+ * \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
+ *
+ */
+
+#include <logog/include/logog.hpp>
+
+// for debugging
+// #include <iostream>
+
+#include "MathLib/LinAlg/BLAS.h"
+
+#include "NonlinearSolver.h"
+
+#include "MathLib/LinAlg/VectorNorms.h"
+
+// TODO: change
+#include "ODETypes.h" // for one shot linear solver
+
+
+namespace NumLib
+{
+
+template<typename Matrix, typename Vector>
+void
+NonlinearSolver<Matrix, Vector, NonlinearSolverTag::Picard>::
+assemble(NonlinearSystem<Matrix, Vector, NonlinearSolverTag::Picard> &sys,
+ Vector const& x) const
+{
+ sys.assembleMatricesPicard(x);
+}
+
+template<typename Matrix, typename Vector>
+bool
+NonlinearSolver<Matrix, Vector, NonlinearSolverTag::Picard>::
+solve(NonlinearSystem<Matrix, Vector, NonlinearSolverTag::Picard> &sys, Vector &x)
+{
+ namespace BLAS = MathLib::BLAS;
+
+ bool success = false;
+
+ for (unsigned iteration=1; iteration<_maxiter; ++iteration)
+ {
+ sys.assembleMatricesPicard(x);
+ sys.getA(_A);
+ sys.getRhs(_rhs);
+ // TODO _x_new might not have been initialized
+ sys.applyKnownComponentsPicard(_A, _rhs, _x_new);
+
+ // std::cout << "A:\n" << Eigen::MatrixXd(A) << "\n";
+ // std::cout << "rhs:\n" << rhs << "\n\n";
+
+ oneShotLinearSolve(_A, _rhs, _x_new);
+
+ BLAS::aypx(x, -1.0, _x_new); // x = _x_new - x
+ auto const error = norm(x);
+ // INFO(" picard iteration %u error: %e", iteration, error);
+
+ x = _x_new;
+
+ if (error < _tol) {
+ success = true;
+ break;
+ }
+
+ if (sys.isLinear()) {
+ // INFO(" picard linear system. not looping");
+ success = true;
+ break;
+ }
+ }
+
+ return success;
+}
+
+
+template<typename Matrix, typename Vector>
+void
+NonlinearSolver<Matrix, Vector, NonlinearSolverTag::Newton>::
+assemble(NonlinearSystem<Matrix, Vector, NonlinearSolverTag::Newton> &sys,
+ Vector const& x) const
+{
+ sys.assembleResidualNewton(x);
+}
+
+template<typename Matrix, typename Vector>
+bool
+NonlinearSolver<Matrix, Vector, NonlinearSolverTag::Newton>::
+solve(NonlinearSystem<Matrix, Vector, NonlinearSolverTag::Newton> &sys, Vector &x)
+{
+ namespace BLAS = MathLib::BLAS;
+
+ bool success = false;
+
+ for (unsigned iteration=1; iteration<_maxiter; ++iteration)
+ {
+ sys.assembleResidualNewton(x);
+ sys.getResidual(x, _res);
+
+ // std::cout << " res:\n" << res << std::endl;
+
+ if (norm(_res) < _tol) {
+ success = true;
+ break;
+ }
+
+ sys.assembleJacobian(x);
+ sys.getJacobian(_J);
+ // TODO _minus_delta_x might not have been initialized
+ sys.applyKnownComponentsNewton(_J, _res, _minus_delta_x);
+
+ // std::cout << " J:\n" << Eigen::MatrixXd(J) << std::endl;
+
+ oneShotLinearSolve(_J, _res, _minus_delta_x);
+
+ // auto const dx_norm = _minus_delta_x.norm();
+ // INFO(" newton iteration %u, norm of delta x: %e", iteration, dx_norm);
+
+ BLAS::axpy(x, -1.0, _minus_delta_x);
+
+ if (sys.isLinear()) {
+ // INFO(" newton linear system. not looping");
+ success = true;
+ break;
+ }
+ }
+
+ return success;
+}
+
+}
diff --git a/NumLib/ODESolver/NonlinearSolver.h b/NumLib/ODESolver/NonlinearSolver.h
new file mode 100644
index 0000000000000000000000000000000000000000..02f1ca009d919020ebebadb204e97f6cd75aaa92
--- /dev/null
+++ b/NumLib/ODESolver/NonlinearSolver.h
@@ -0,0 +1,149 @@
+/**
+ * \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 NUMLIB_NONLINEARSOLVER_H
+#define NUMLIB_NONLINEARSOLVER_H
+
+#include <logog/include/logog.hpp>
+
+#include "Types.h"
+#include "NonlinearSystem.h"
+
+
+namespace NumLib
+{
+
+//! \addtogroup ODESolver
+//! @{
+
+/*! Find a solution to a nonlinear equation.
+ *
+ * \tparam Matrix the type of matrices occuring in the linearization of the equation.
+ * \tparam Vector the type of the solution vector of the equation.
+ * \tparam NLTag a tag indicating the method used for solving the equation.
+ */
+template<typename Matrix, typename Vector, NonlinearSolverTag NLTag>
+class NonlinearSolver;
+
+
+/*! Find a solution to a nonlinear equation using the Newton-Raphson method.
+ *
+ * \tparam Matrix the type of matrices occuring in the linearization of the equation.
+ * \tparam Vector the type of the solution vector of the equation.
+ */
+template<typename Matrix, typename Vector>
+class NonlinearSolver<Matrix, Vector, NonlinearSolverTag::Newton> final
+{
+public:
+ //! Type of the nonlinear equation system to be solved.
+ using System = NonlinearSystem<Matrix, Vector, NonlinearSolverTag::Newton>;
+
+ /*! Constructs a new instance.
+ *
+ * \param tol the tolerance of the solver. \todo Be more specific about that!
+ * \param maxiter the maximum number of iterations used to solve the equation.
+ */
+ explicit
+ NonlinearSolver(double const tol, const unsigned maxiter)
+ : _tol(tol)
+ , _maxiter(maxiter)
+ {}
+
+ /*! Only assemble the equation system.
+ *
+ * \note This method is needed to preload CrankNicolson time discretization scheme.
+ * It is not used for the general solver steps; in those only the solve() method
+ * is needed.
+ *
+ * \param sys the equation system to be assembled.
+ * \param x the state at which the equation system will be assembled.
+ */
+ void assemble(System& sys, Vector const& x) const;
+
+ /*! Assemble and solve the equation system.
+ *
+ * \param sys the equation system to be solved.
+ * \param x in: the initial guess, out: the solution.
+ *
+ * \retval true if the equation system could be solved
+ * \retval false otherwise
+ */
+ bool solve(System& sys, Vector& x);
+
+private:
+ const double _tol; //!< tolerance of the solver
+ const unsigned _maxiter; //!< maximum number of iterations
+
+ Vector _res; //!< The residual.
+ Matrix _J; //!< The Jacobian of the residual.
+ Vector _minus_delta_x; //!< The Newton-Raphson method solves the linearized equation
+ //!< \f$ J \cdot (-\Delta x) = r \f$ repeatedly.
+};
+
+
+/*! Find a solution to a nonlinear equation using the Picard fixpoint iteration method.
+ *
+ * \tparam Matrix the type of matrices occuring in the linearization of the equation.
+ * \tparam Vector the type of the solution vector of the equation.
+ */
+template<typename Matrix, typename Vector>
+class NonlinearSolver<Matrix, Vector, NonlinearSolverTag::Picard> final
+{
+public:
+ //! Type of the nonlinear equation system to be solved.
+ using System = NonlinearSystem<Matrix, Vector, NonlinearSolverTag::Picard>;
+
+ /*! Constructs a new instance.
+ *
+ * \param tol the tolerance of the solver. \todo Be more specific about that!
+ * \param maxiter the maximum number of iterations used to solve the equation.
+ */
+ explicit
+ NonlinearSolver(double const tol, const unsigned maxiter)
+ : _tol(tol)
+ , _maxiter(maxiter)
+ {}
+
+ /*! Only assemble the equation system.
+ *
+ * \note This method is needed to preload CrankNicolson time discretization scheme.
+ * It is not used for the general solver steps; in those only the solve() method
+ * is needed.
+ *
+ * \param sys the equation system to be assembled.
+ * \param x the state at which the equation system will be assembled.
+ */
+ void assemble(System& sys, Vector const& x) const;
+
+ /*! Assemble and solve the equation system.
+ *
+ * \param sys the equation system to be solved.
+ * \param x in: the initial guess, out: the solution.
+ *
+ * \retval true if the equation system could be solved
+ * \retval false otherwise
+ */
+ bool solve(System& sys, Vector& x);
+
+private:
+ const double _tol; //!< tolerance of the solver
+ const unsigned _maxiter; //!< maximum number of iterations
+
+ Matrix _A; //!< \c Matrix describing the linearized system.
+ Vector _rhs; //!< \c Vector describing the linearized system.
+ Vector _x_new; //!< \c Vector to store solutions of \f$ A \cdot x = \mathit{rhs} \f$.
+};
+
+//! @}
+
+}
+
+#include "NonlinearSolver-impl.h"
+
+#endif // NUMLIB_NONLINEARSOLVER_H