--- a/+rv/+time/RungekuttaExteriorRvBdf.m	Tue Mar 05 10:57:26 2019 +0100
+++ b/+rv/+time/RungekuttaExteriorRvBdf.m	Wed Mar 06 09:45:52 2019 +0100
@@ -5,7 +5,8 @@
         t       % Time point
         v       % Solution vector
         n       % Time level
-        coeffs  % The coefficents used for the RK time integration
+        rkScheme  % The particular RK scheme used for time integration
+
         
         % Properties related to the residual viscositys
         RV              % Residual Viscosity operator
@@ -15,6 +16,8 @@
                         % dictates which accuracy the boot-strapping should start from.
         upperBdfOrder   % Orders of the approximation of the time deriative, used for the RV evaluation.
                         % Dictates the order of accuracy used once the boot-strapping is complete.
+
+
     end
     methods
         function obj = RungekuttaExteriorRvBdf(F, k, t0, v0, RV, rkOrder, bdfOrders)
@@ -23,17 +26,23 @@
             obj.t = t0;
             obj.v = v0;
             obj.n = 0;
-            % Extract the coefficients for the specified rkOrder
-            % used for the RK updates from the Butcher tableua.
-            [s,a,b,c] = time.rk.butcherTableau(rkOrder);
-            obj.coeffs = struct('s',s,'a',a,'b',b,'c',c);
-        
             obj.RV = RV;
             obj.lowerBdfOrder = bdfOrders.lowerBdfOrder;
             obj.upperBdfOrder = bdfOrders.upperBdfOrder;
             assert((obj.lowerBdfOrder >= 1) && (obj.upperBdfOrder <= 6));
             obj.v_prev = [];
             obj.DvDt = rv.time.BdfDerivative();
+
+            if (rkOrder == 4) % Use specialized RK4 scheme
+                obj.rkScheme = @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(rkOrder);
+                coeffs = struct('s',s,'a',a,'b',b,'c',c);
+                obj.rkScheme = @(v,t,dt,F) time.rk.rungekutta(v, t , dt, F, coeffs);
+            end
+
         end
 
         function [v, t] = getV(obj)
@@ -74,8 +83,9 @@
             end
 
             % Fix the viscosity of the RHS function F
-            F_visc = @(v,t) obj.F(v, t, viscosity);
-            obj.v = time.rk.rungekutta(obj.v, obj.t, obj.k, F_visc, obj.coeffs);
+            m = length(viscosity);
+            F_visc = @(v,t) obj.F(v, t, spdiags(viscosity,0,m,m));
+            obj.v = obj.rkScheme(obj.v, obj.t, obj.k, F_visc);
             obj.t = obj.t + obj.k;
             obj.n = obj.n + 1;
         end