--- a/+time/Rungekutta.m	Thu Sep 20 19:02:01 2018 +0200
+++ b/+time/Rungekutta.m	Fri Sep 21 09:13:03 2018 +0200
@@ -5,7 +5,7 @@
         t       % Time point
         v       % Solution vector
         n       % Time level
-        coeffs  % The coefficents used for the RK time integration
+        scheme  % The coefficents used for the RK time integration
     end
 
 
@@ -19,10 +19,17 @@
             obj.t = t0;
             obj.v = v0;
             obj.n = 0;
-            % Extract the coefficients for the specified order
-            % used for the RK updates from the Butcher tableua.
-            [s,a,b,c] = butcherTableau(order);
-            obj.coeffs = struct('s',s,'a',a,'b',b,'c',c);
+            % TBD: Order 4 is also implemented in the butcher tableau, but the rungekutta_4.m implementation
+            % might be slightly more efficient. Need to do some profiling before deciding whether or not to keep it.
+            if (order == 4)
+                obj.scheme = @time.rk.rungekutta_4;
+            else
+                % Extract the coefficients for the specified order
+                % used for the RK updates from the Butcher tableua.
+                [s,a,b,c] = time.rk.butcherTableau(order);
+                coeffs = struct('s',s,'a',a,'b',b,'c',c);
+                obj.scheme = @(v,t,dt,F) time.rk.rungekutta(v, t , dt, F, coeffs);
+            end
         end
 
         function [v,t] = getV(obj)
@@ -31,7 +38,7 @@
         end
 
         function obj = step(obj)
-            obj.v = time.rk.rungekutta(obj.v, obj.t, obj.k, obj.F, obj.coeffs);
+            obj.v = obj.scheme(obj.v, obj.t, obj.k, obj.F);
             obj.t = obj.t + obj.k;
             obj.n = obj.n + 1;
         end