--- a/+rv/+time/RungekuttaExteriorRv.m	Tue Mar 05 10:57:26 2019 +0100
+++ b/+rv/+time/RungekuttaExteriorRv.m	Wed Mar 06 09:45:52 2019 +0100
@@ -5,7 +5,7 @@
         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
         RV              % Residual Viscosity operator
         DvDt            % Function for computing the time deriative used for the RV evaluation
     end
@@ -17,10 +17,16 @@
             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] = time.rk.butcherTableau(order);
-            obj.coeffs = struct('s',s,'a',a,'b',b,'c',c);
+            
+            if (order == 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(order);
+                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
         
             obj.RV = RV;
             obj.DvDt = DvDt;
@@ -41,8 +47,10 @@
         % obj.coeffs, using a fixed residual viscosity for the Runge-Kutta substeps
         function obj = step(obj)            
             % Fix the viscosity of the RHS function F
-            F_visc = @(v,t) obj.F(v,t,obj.RV.evaluateViscosity(obj.v, obj.DvDt(obj.v)));
-            obj.v = time.rk.rungekutta(obj.v, obj.t, obj.k, F_visc, obj.coeffs);
+            viscosity = obj.RV.evaluateViscosity(obj.v, obj.DvDt(obj.v));
+            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