diff --git a/src/Makefile.am b/src/Makefile.am
index 4adfdcb1ae423f7a607478317ed25aa36ba32ffc..5d7f6d03738a5046e0246d0ec9325024d1b1ac79 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -5,6 +5,12 @@ check_PROGRAMS = \
 	test-python \
 	test-gradient-method
 
+bin_PROGRAMS = \
+	gauss-seidel-sample
+
+gauss_seidel_sample_SOURCES = \
+	gauss-seidel-sample.cc gaussseidel.hh
+
 test_python_SOURCES = \
 	test-python.cc
 
diff --git a/src/gauss-seidel-sample.cc b/src/gauss-seidel-sample.cc
new file mode 100644
index 0000000000000000000000000000000000000000..cd60bc7624ddbde8155461ce13706231b75a9f26
--- /dev/null
+++ b/src/gauss-seidel-sample.cc
@@ -0,0 +1,50 @@
+/* -*- mode:c++; mode:semantic -*- */
+
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <dune/common/fmatrix.hh>
+#include <dune/common/fvector.hh>
+
+#include "nicefunction.hh"
+#include "mynonlinearity.hh"
+
+#include "gaussseidel.hh"
+
+int main() {
+  size_t const dimension = 3;
+  typedef Dune::FieldMatrix<double, dimension, dimension> LocalMatrix;
+  typedef Dune::FieldVector<double, dimension> LocalVector;
+  size_t const nodes = 10;
+  typedef Dune::FieldMatrix<LocalMatrix, nodes, nodes> GlobalMatrix;
+  typedef Dune::FieldVector<LocalVector, nodes> GlobalVector;
+
+  GlobalMatrix A;
+  // semi-random data
+  for (size_t i = 0; i < nodes; ++i)
+    for (size_t k = 0; k < nodes; ++k)
+      for (size_t j = 0; j < dimension; ++j)
+        for (size_t l = 0; l < dimension; ++l)
+          A[i][k][j][l] = ((i == k) && (j == l) ? 1 : 0);
+
+  GlobalVector b;
+  // semi-random data
+  for (size_t i = 0; i < nodes; ++i)
+    for (size_t j = 0; j < dimension; ++j)
+      b[i][j] = 1;
+
+  Dune::SampleFunction const f = Dune::SampleFunction(); // semi-random function
+  Dune::MyNonlinearity<dimension> const phi(f);
+
+  GlobalVector unu(LocalVector(0));
+  // semi-random data
+  for (size_t m = 0; m < nodes; ++m)
+    for (size_t j = 0; j < dimension; ++j)
+      unu[m][j] = m + j;
+
+  GaussSeidel(unu, A, b, phi);
+  std::cout << std::endl << std::endl;
+  GaussSeidel(unu, A, b, phi);
+  return 0;
+}
diff --git a/src/gaussseidel.hh b/src/gaussseidel.hh
new file mode 100644
index 0000000000000000000000000000000000000000..46db7f1cff5900553f98c71a724e5f42ee479101
--- /dev/null
+++ b/src/gaussseidel.hh
@@ -0,0 +1,52 @@
+/* -*- mode:c++; mode:semantic -*- */
+
+#include <dune/common/fmatrix.hh>
+#include <dune/common/fvector.hh>
+
+#include "nicefunction.hh"
+#include "mynonlinearity.hh"
+#include "samplefunctional.hh"
+
+// Just for debugging
+template <int dim, int nodes>
+double J_eval(
+    Dune::FieldMatrix<Dune::FieldMatrix<double, dim, dim>, nodes, nodes> const &
+        A,
+    Dune::FieldVector<Dune::FieldVector<double, dim>, nodes> const &b,
+    Dune::MyNonlinearity<dim> phi,
+    Dune::FieldVector<Dune::FieldVector<double, dim>, nodes> const &v) {
+  double ret = 0;
+  for (size_t i = 0; i < nodes; ++i) {
+    for (size_t k = 0; k < nodes; ++k) {
+      Dune::FieldMatrix<double, dim, dim> const localA = A[i][k];
+      Dune::FieldVector<double, dim> tmp;
+      localA.mv(v[k], tmp);
+      ret += 1 / 2 * (tmp * v[i]);
+    }
+    ret -= b[i] * v[i];
+    ret +=
+        phi(v[i]); // an example of phi^i(v[i]) (phi is independent of i here)
+  }
+  return ret;
+}
+
+template <int dim, int nodes>
+void GaussSeidel(
+    Dune::FieldVector<Dune::FieldVector<double, dim>, nodes> &unu,
+    Dune::FieldMatrix<Dune::FieldMatrix<double, dim, dim>, nodes, nodes> const &
+        A,
+    Dune::FieldVector<Dune::FieldVector<double, dim>, nodes> const &b,
+    Dune::MyNonlinearity<dim> const &phi) {
+  for (size_t m = 0; m < nodes; ++m) {
+    typename Dune::SampleFunctional<dim>::SmallVector const localb = b[m];
+    typename Dune::SampleFunctional<dim>::SmallMatrix const localA = A[m][m];
+    Dune::SampleFunctional<dim> J(localA, localb, phi);
+    typename Dune::SampleFunctional<dim>::SmallVector correction;
+    for (size_t local_iteration = 0; local_iteration < 2;
+         ++local_iteration) { // 2 is random here
+      Dune::minimise(J, unu[m], correction);
+      unu[m] += correction;
+    }
+    std::cout << J_eval(A, b, phi, unu) << std::endl; // debugging
+  }
+}