diff --git a/python/dune/codegen/options.py b/python/dune/codegen/options.py
index 72697492907cbf4afdfe62d1a3789606a8e2c290..51722da3fe34a3fa24b85f4d47e06750c0cc0007 100644
--- a/python/dune/codegen/options.py
+++ b/python/dune/codegen/options.py
@@ -37,6 +37,7 @@ class CodegenGlobalOptionsArray(ImmutableRecord):
     debug_cache_with_stack = CodegenOption(default=False, helpstr="Store stack along with cache objects. Makes debugging caching issues easier.")
     driver_file = CodegenOption(helpstr="The filename for the generated driver header")
     explicit_time_stepping = CodegenOption(default=False, helpstr="use explicit time stepping")
+    time_stepping_order = CodegenOption(default=1, helpstr="Order of the time stepping method")
     exact_solution_expression = CodegenOption(helpstr="name of the exact solution expression in the ufl file")
     compare_l2errorsquared = CodegenOption(helpstr="maximal allowed l2 error squared of difference between numerical solution and interpolation of exact solution (NOTE: requires --exact-solution-expression)")
     grid_info = CodegenOption(default=None, helpstr="Path to file with information about facedir and facemod variations. This can be used to limit the generation of skeleton kernels.")
diff --git a/python/dune/codegen/pdelab/driver/instationary.py b/python/dune/codegen/pdelab/driver/instationary.py
index 5cf7170f1dca6f4274db4bf77b4b1627466e20cb..355fd0e743f6a59716319241018d3c5794c784bd 100644
--- a/python/dune/codegen/pdelab/driver/instationary.py
+++ b/python/dune/codegen/pdelab/driver/instationary.py
@@ -132,10 +132,25 @@ def name_time():
 def typedef_timesteppingmethod(name):
     r_type = type_range()
     explicit = get_option('explicit_time_stepping')
+    order = get_option('time_stepping_order')
     if explicit:
-        return "using {} = Dune::PDELab::ExplicitEulerParameter<{}>;".format(name, r_type)
+        if order == 1:
+            return "using {} = Dune::PDELab::ExplicitEulerParameter<{}>;".format(name, r_type)
+        elif order == 2:
+            return "using {} = Dune::PDELab::HeunParameter<{}>;".format(name, r_type)
+        elif order == 3:
+            return "using {} = Dune::PDELab::Shu3Parameter<{}>;".format(name, r_type)
+        elif order == 4:
+            return "using {} = Dune::PDELab::RK4Parameter<{}>;".format(name, r_type)
+        else:
+            raise NotImplementedError("Time stepping method not supported")
     else:
-        return "using {} = Dune::PDELab::OneStepThetaParameter<{}>;".format(name, r_type)
+        if order == 1:
+            return "using {} = Dune::PDELab::OneStepThetaParameter<{}>;".format(name, r_type)
+        elif order == 2:
+            return "using {} = Dune::PDELab::Alexander2Parameter<{}>;".format(name, r_type)
+        elif order == 3:
+            return "using {} = Dune::PDELab::Alexander3Parameter<{}>;".format(name, r_type)
 
 
 def type_timesteppingmethod():
@@ -150,8 +165,12 @@ def define_timesteppingmethod(name):
     if explicit:
         return "{} {};".format(tsm_type, name)
     else:
-        ini = name_initree()
-        return "{} {}({}.get<double>(\"instat.theta\",1.0));".format(tsm_type, name, ini)
+        order = get_option('time_stepping_order')
+        if order == 1:
+            ini = name_initree()
+            return "{} {}({}.get<double>(\"instat.theta\",1.0));".format(tsm_type, name, ini)
+        else:
+            return "{} {};".format(tsm_type, name)
 
 
 def name_timesteppingmethod():