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