diff --git a/extensions/demo/other/src/polarizationField.cc b/extensions/demo/other/src/polarizationField.cc
new file mode 100644
index 0000000000000000000000000000000000000000..46af42e40f1dac44a5d6a51d53a5034635640a38
--- /dev/null
+++ b/extensions/demo/other/src/polarizationField.cc
@@ -0,0 +1,99 @@
+/******************************************************************************
+ *
+ * Extension of AMDiS - Adaptive multidimensional simulations
+ *
+ * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
+ * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
+ *
+ * Authors: Simon Praetorius et al.
+ *
+ * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
+ * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+ *
+ *
+ * See also license.opensource.txt in the distribution.
+ * 
+ ******************************************************************************/
+
+#include "AMDiS.h"
+#include "PolarizationField_RB.h"
+#include "time/ExtendedRosenbrockAdaptInstationary.h"
+   
+using namespace AMDiS;
+
+class MyPolarizationField : public PolarizationField_RB
+{
+public:
+  MyPolarizationField(std::string name)
+  : PolarizationField_RB(name) {}
+  
+  /// Set initial condition and perform initial refinement
+  void solveInitialProblem(AdaptInfo *adaptInfo)
+  {
+    struct RandomNormalizedVector : public AbstractFunction<WorldVector<double>, WorldVector<double> >
+    {
+      RandomNormalizedVector() : AbstractFunction<WorldVector<double>, WorldVector<double> >(1) 
+      {
+	std::srand(time(0));
+      }
+      
+      WorldVector<double> operator()(const WorldVector<double>& x) const
+      {
+	WorldVector<double> p;
+	for (int i = 0; i < p.getSize(); i++)
+	  p[i] = cos(4.0*m_pi * ((std::rand() / static_cast<double>(RAND_MAX)) - 0.5));
+	
+	double nrm = norm(p);
+	for (int i = 0; i < p.getSize(); i++)
+	  p[i] *= 1.0/nrm;
+	return p;
+      }
+    };
+    
+    self::getVectorField()->interpol(new RandomNormalizedVector);
+    for (int i = 0; i < Global::getGeo(WORLD); i++)
+      transformDOF(self::getVectorField(), self::getSolution()->getDOFVector(i), new AMDiS::Component2<>(i));
+  }
+  
+  void fillBoundaryConditions() override
+  {
+    std::vector<DOFVector<double>*> values(3);
+    values[0] = new DOFVector<double>(self::getFeSpace(0), "v0"); values[0]->set(0.0);
+    values[1] = new DOFVector<double>(self::getFeSpace(0), "v1"); values[1]->set(1.0);
+    values[2] = new DOFVector<double>(self::getFeSpace(0), "v2"); values[2]->set(-1.0);
+    int idx = 1;
+    for (int i = 0; i < Global::getGeo(WORLD); i++) {
+      for (int j = 0; j < Global::getGeo(WORLD); j++) {
+	double value1 = double(i == j);
+	double value2 = -double(i == j);
+	self::prob->addDirichletBC(idx, 1-j + self::dow, 1-j, new AMDiS::Const<double, WorldVector<double> >(value1));
+	self::prob->addDirichletBC(idx+1, 1-j + self::dow, 1-j, new AMDiS::Const<double, WorldVector<double> >(value1));
+      }
+      idx += 2;
+    }
+  }
+};
+
+
+
+int main(int argc, char** argv)
+{ FUNCNAME("main");
+
+  AMDiS::init(argc, argv);  
+  Timer t;
+  
+  MyPolarizationField pProb("p");
+  pProb.initialize(INIT_ALL);
+  pProb.initTimeInterface();
+
+  // Adapt-Infos
+  AdaptInfo adaptInfo("adapt", pProb.getNumComponents());
+  ExtendedRosenbrockAdaptInstationary<MyPolarizationField> adaptInstat("adapt", pProb, adaptInfo, pProb, adaptInfo);
+
+  int error_code = adaptInstat.adapt(); 
+
+  MSG("elapsed time= %d sec\n", t.elapsed());
+  AMDiS::finalize();
+
+  return error_code;
+}