--- 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

--- 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

--- a/+rv/constructDiffOps.m	Tue Mar 05 10:57:26 2019 +0100
+++ b/+rv/constructDiffOps.m	Wed Mar 06 09:45:52 2019 +0100
@@ -1,11 +1,11 @@
 function [D_rv, D_flux, DvDt, solutionPenalties, residualPenalties] = constructDiffOps(scheme, g, order, schemeParams, opSet, BCs)
     %% DiffOps for solution vector
-    [D, solutionPenalties] = constructTotalFluxDiffOp(scheme, g, order, schemeParams, opSet, BCs);
+    [D, solutionPenalties] = constructFluxDiffOp(scheme, g, order, schemeParams, opSet, BCs);
     D2 = constructSymmetricD2Operator(g, order, opSet);
-    D_rv = @(v,viscosity)(D(v) + D2(v, viscosity));
+    D_rv = @(v,Viscosity)(D(v) + D2(Viscosity)*v);
 
     %% DiffOps for residual viscosity
-    [D_flux, residualPenalties] = constructTotalFluxDiffOp(scheme, g, max(order-2,2), schemeParams, opSet, BCs);
+    [D_flux, residualPenalties] = constructFluxDiffOp(scheme, g, max(order-2,2), schemeParams, opSet, BCs);
     % DiffOp for flux in residual viscosity. Due to sign conventions of the implemnted schemes, we need to
     % change the sign.
     D_flux = @(v) -D_flux(v);
@@ -13,7 +13,7 @@
     DvDt = D;
 end
 
-function [D, penalties] = constructTotalFluxDiffOp(scheme, g, order, schemeParams, opSet, BCs)
+function [D, penalties] = constructFluxDiffOp(scheme, g, order, schemeParams, opSet, BCs)
     diffOp = scheme(g, order, schemeParams{:}, opSet);
     [D, penalties]  = addClosuresToDiffOp(diffOp, BCs);
 end
@@ -48,10 +48,8 @@
        I{i} = speye(m(i));
     end
 
-    % TBD: How is this generalized to a loop over dimensions or similar?
     switch g.D()
         case 1
-            
             e_r = ops{1}.e_r;
             e_l = ops{1}.e_l;
             Hi = ops{1}.HI;
@@ -59,9 +57,11 @@
             if isequal(opSet,@sbp.D1Upwind)
                 Dm = ops{1}.Dm;
                 Dp = ops{1}.Dp;
-                D2 = @(viscosity) Dm*spdiag(viscosity)*Dp-Hi*(B*spdiag(viscosity)*Dp);
+                M =  Dm - Hi*B;
+                D2 = @(Viscosity) M*Viscosity*Dp;
             else
-                D2 = @(viscosity)ops{1}.D2(viscosity);
+                % TODO: Fix closure for D2Variable
+                D2 = @(Viscosity)ops{1}.D2(Viscosity);
             end
         case 2
             % TODO: 
@@ -75,7 +75,8 @@
             Dp_x  = kron(ops{1}.Dp,I{2});
             H_x  = kron(ops{1}.HI,I{2});
             B_x = e_e*e_e' - e_w*e_w';
-            D2_x = @(viscosity) Dm_x*spdiag(viscosity)*Dp_x-H_x*(B_x*spdiag(viscosity)*Dp_x);
+            M_x = Dm_x-H_x*B_x;
+            D2_x = @(Viscosity) M_x*Viscosity*Dp_x;
 
             e_n = kron(I{1},ops{2}.e_r);
             e_s = kron(I{1},ops{2}.e_l);
@@ -83,10 +84,10 @@
             Dp_y  = kron(I{1},ops{2}.Dp);
             H_y = kron(I{1},ops{2}.HI);
             B_y = e_n*e_n' - e_s*e_s';
-            D2_y = @(viscosity) Dm_y*spdiag(viscosity)*Dp_y-H_y*(B_y*spdiag(viscosity)*Dp_y);
-            D2 = @(viscosity)D2_x(viscosity) + D2_y(viscosity);
+            M_y = Dm_y-H_y*B_y;
+            D2_y = @(Viscosity) M_y*Viscosity*Dp_y;
+            D2 = @(Viscosity) D2_x(Viscosity) + D2_y(Viscosity);
         otherwise
             error('3D not yet implemented')
     end
-    D2 = @(v, viscosity) D2(viscosity)*v;
 end
\ No newline at end of file

--- a/+scheme/Burgers2d.m	Tue Mar 05 10:57:26 2019 +0100
+++ b/+scheme/Burgers2d.m	Wed Mar 06 09:45:52 2019 +0100
@@ -47,25 +47,28 @@
             if (isequal(opSet,@sbp.D1Upwind))
                 Dx = kron((ops_x.Dp + ops_x.Dm)/2,Iy);
                 Dy = kron(Ix,(ops_y.Dp + ops_y.Dm)/2);
-                DissOpx = kron((ops_x.Dm - ops_x.Dp)/2,Iy);
-                DissOpy = kron(Ix,(ops_y.Dm - ops_y.Dp)/2);
+                DissOpx = kron((ops_x.Dp - ops_x.Dm)/2,Iy);
+                DissOpy = kron(Ix,(ops_y.Dp - ops_y.Dm)/2);
                 D1 = Dx + Dy;   
                 switch pde_form
                     case 'skew-symmetric'
+                        D = -1/3*D1
                         switch length(fluxSplitting)
                             case 1
                                 DissOp = DissOpx + DissOpy;
-                                obj.D = @(v) -(1/3*D1*v.*v + (1/3*spdiag(v)*D1 + fluxSplitting{1}(v)*DissOp)*v);
+                                obj.D = @(v) D*(v.*v) + (spdiags(v,0,m_tot,m_tot)*D + fluxSplitting{1}(v)*DissOp)*v;
                             case 2
-                                obj.D = @(v) -(1/3*D1*v.*v + (1/3*spdiag(v)*D1 + fluxSplitting{1}(v)*DissOpx + fluxSplitting{2}(v)*DissOpy)*v);
+                                obj.D = @(v) D*(v.*v) + (spdiags(v,0,m_tot,m_tot)*D + fluxSplitting{1}(v)*DissOpx + fluxSplitting{2}(v)*DissOpy)*v;
                         end
                     case 'conservative'
+                        D = -1/2*D1;
                         switch length(fluxSplitting)
                             case 1
                                 DissOp = DissOpx + DissOpy;
-                                obj.D = @(v) -(1/2*D1*v.*v + fluxSplitting{1}(v)*DissOp*v);
+                                % TODO: Check if we can use fluxSplitting{1} here instead
+                                obj.D = @(v) D*(v.*v) + max(abs(v))*DissOp*v;
                             case 2
-                                obj.D = @(v) -(1/2*D1*v.*v + (fluxSplitting{1}(v)*DissOpx + fluxSplitting{2}(v)*DissOpy)*v);
+                                obj.D = @(v) D*(v.*v) + (fluxSplitting{1}(v)*DissOpx + fluxSplitting{2}(v)*DissOpy)*v;
                         end
                         
                 end
@@ -75,9 +78,11 @@
                 D1 = Dx + Dy;
                 switch pde_form
                     case 'skew-symmetric'
-                        obj.D = @(v) -(1/3*D1*v.*v + 1/3*spdiag(v)*D1*v);
+                        D = -1/3*D1
+                        obj.D = @(v) D*(v.*v) + spdiags(v,0,m_tot,m_tot)*D*v;
                     case 'conservative'
-                        obj.D = @(v) -1/2*D1*v.*v;
+                        D = -1/2*D1;
+                        obj.D = @(v) D*(v.*v);
                 end
             end
 
@@ -103,14 +108,11 @@
                 case {'D', 'd', 'dirichlet', 'Dirichlet'}
 
                     magnitude = 2/3;
-                    tau = @(v) s*magnitude*obj.Hi*e*H_b*spdiag((v(index)-s*abs(v(index)))/2);
-                    closure = @(v) tau(v)*v(index);
+                    Tau = s*magnitude*obj.Hi*e*H_b/2;
+                    m = length(index);
+                    tau = @(v) Tau*spdiags((v(index)-s*abs(v(index))),0,m,m);
+                    closure = @(v) Tau*((v(index)-s*abs(v(index))).*v(index));
                     penalty = @(v) -tau(v);
-
-
-                    % tau = s*e*H_b;
-                    % closure = @(v) (obj.Hi*tau)*(((v(index)-s*abs(v(index)))/3).*v(index));
-                    % penalty = -obj.Hi*tau;
                 otherwise
                     error('No such boundary condition: type = %s',type);
             end