diff -r e3d8f24b2c1c -r cc2fd58989d1 +rv/+time/getRvTimestepper.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+rv/+time/getRvTimestepper.m	Fri Jun 28 17:24:01 2019 +0200
@@ -0,0 +1,57 @@
+function ts = getRvTimestepper(opt, diffOpStruct, residualViscosity, data, t0, v0)
+    switch opt.method
+        case 'rkRvBdf'
+            ts = rkRvBdfTimestepper(opt, diffOpStruct, residualViscosity, data, t0, v0);
+        case 'rkRvMs'
+            ts = rkRvMsTimestepper(opt, diffOpStruct, residualViscosity, data, t0, v0);
+        case 'rkRvMg'
+            ts = rkRvMgTimestepper(opt, diffOpStruct, residualViscosity, data, t0, v0);
+        case 'rkRvInstage'
+            ts = rkRvInstageTimestepper(opt, diffOpStruct, residualViscosity, data, t0, v0);
+        otherwise
+            error('Timestepping method ''%s'' not supported',method);
+    end
+end
+
+function F = stabilizedRhs(D, data)
+    if isa(data, 'function_handle')
+        F = @(v, t, viscosity) D(v, viscosity) + data(v);
+    else
+        F = @(v, t, viscosity) D(v, viscosity) + data;
+    end
+end
+
+function F = unstabilizedRhs(D, data)
+    if isa(data, 'function_handle')
+        F = @(v, t) D(v) + data(v);
+    else
+        F = @(v, t) D(v) + data;
+    end
+end
+
+function ts = rkRvBdfTimestepper(opt, diffOpStruct, residualViscosity, data, t0, v0)
+    F = stabilizedRhs(diffOpStruct.D_scheme, data);
+    ts = rv.time.RungekuttaRvBdf(F, opt.k, t0, v0, residualViscosity, opt.rkOrder, opt.bdfOrders);
+end
+
+function ts = rkRvMsTimestepper(opt, diffOpStruct, residualViscosity, data, t0, v0)
+    F = stabilizedRhs(diffOpStruct.D_scheme, data);
+    F_unstab = unstabilizedRhs(diffOpStruct.D_unstable, data);
+    ts = rv.time.RungekuttaRvMultiStage(F, F_unstab, opt.k, t0, v0,...
+                                        residualViscosity, diffOpStruct.D_t, opt.rkOrder);
+end
+
+function ts = rkRvMgTimestepper(opt, diffOpStruct, residualViscosity, data, t0, v0)
+    F = stabilizedRhs(diffOpStruct.D_scheme, data);
+    F_coarse = unstabilizedRhs(diffOpStruct.D_coarse, data);
+    ts = rv.time.RungekuttaRvMultiGrid(F, F_coarse, opt.k, t0, v0,...
+                                       residualViscosity, diffOpStruct.D_t, opt.rkOrder);
+end
+
+function ts = rkRvInstageTimestepper(opt, diffOpStruct, residualViscosity, data, t0, v0)
+    F = stabilizedRhs(diffOpStruct.D_scheme, data);
+    ts = rv.time.RungekuttaRvInstage(F, opt.k, t0, v0, residualViscosity, diffOpStruct.D_t, opt.rkOrder);
+end
+
+
+