diff -r fbb8be3177c8 -r cda996e64925 +sbp/dissipationOperator.m
--- a/+sbp/dissipationOperator.m	Thu Sep 27 09:30:21 2018 +0200
+++ b/+sbp/dissipationOperator.m	Fri Oct 12 08:41:57 2018 +0200
@@ -63,6 +63,5 @@
         otherwise
             error('Order not yet supported', order);
     end
-    scaling
-    D = scaling*Hinv*D;
+    D = scaling*sparse(Hinv*D);
 end

diff -r fbb8be3177c8 -r cda996e64925 +scheme/Burgers1D.m
--- a/+scheme/Burgers1D.m	Thu Sep 27 09:30:21 2018 +0200
+++ b/+scheme/Burgers1D.m	Fri Oct 12 08:41:57 2018 +0200
@@ -20,10 +20,6 @@
             assert(grid.size() == length(params.eps));
             m = grid.size();
             lim = grid.lim{1}; % Ugly, and only applicable for cartesian grids.
-            default_arg('pde_form','skew-symmetric');
-            default_arg('operator_type','narrow');
-            default_arg('dissipation','on');
-
             switch operator_type
                 case 'narrow'
                     ops = sbp.D4Variable(m, lim, order);
@@ -65,7 +61,6 @@
                     error('Other operator types not yet supported', operator_type);
             end
 
-            %% TODO: Figure out how to evaluate viscosity as viscosity(v,t) here instead of parametrizing D on the viscosity.
             switch pde_form
                 case 'skew-symmetric'
                     if (strcmp(dissipation,'on'))
@@ -79,6 +74,8 @@
                     else
                         D = @(v, viscosity) -1/2*D1*v.^2 + D2(params.eps + viscosity)*v;
                     end
+                otherwise
+                    error('Not supported', pde_form);
             end
 
             obj.grid = grid;
@@ -94,17 +91,15 @@
             obj.d_r =  d_r;
         end
 
-        % Closure functions return the opertors applied to the own doamin to close the boundary
-        % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other domain.
+        % Closure functions return the operators applied to the own doamin to close the boundary
+        % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other domain.
         %       boundary            is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
         %       type                is a string specifying the type of boundary condition if there are several.
         %       data                is a function returning the data that should be applied at the boundary.
-        %       neighbour_scheme    is an instance of Scheme that should be interfaced to.
-        %       neighbour_boundary  is a string specifying which boundary to interface to.
         function [closure, penalty] = boundary_condition(obj,boundary,type,data)
             default_arg('type','robin');
             default_arg('data',0);
-            [e, s, d, i_b] = obj.get_boundary_ops(boundary);
+            [e, d, i_b, s] = obj.get_boundary_ops(boundary);
             switch type
                 % Stable robin-like boundary conditions ((u+-abs(u))*u/3 - eps*u_x)) with +- at left/right boundary
                 case {'R','robin'}
@@ -118,26 +113,25 @@
                         otherwise
                             error('Wierd data argument!')
                     end
-                % Unknown, boundary condition
                 otherwise
                     error('No such boundary condition: type = %s',type);
             end
         end
 
-        % Ruturns the boundary ops and sign for the boundary specified by the string boundary.
+        % Ruturns the boundary ops, boundary index and sign for the boundary specified by the string boundary.
         % The right boundary is considered the positive boundary
-        function [e, s, d, i_b] = get_boundary_ops(obj,boundary)
+        function [e, d, i_b, s] = get_boundary_ops(obj,boundary)
             switch boundary
                 case 'l'
                     e = obj.e_l;
-                    s = -1;
                     d = obj.d_l;
                     i_b = 1;
+                    s = -1;
                 case 'r'
                     e = obj.e_r;
-                    s = 1;
                     d = obj.d_r;
                     i_b = length(e);
+                    s = 1;
                 otherwise
                     error('No such boundary: boundary = %s',boundary);
             end

diff -r fbb8be3177c8 -r cda996e64925 +time/RungekuttaRV.m
--- a/+time/RungekuttaRV.m	Thu Sep 27 09:30:21 2018 +0200
+++ b/+time/RungekuttaRV.m	Fri Oct 12 08:41:57 2018 +0200
@@ -29,17 +29,14 @@
         end
 
         function state = getState(obj)
-            [residual, u_t, f_x] = obj.RV.getResidual();
-            state = struct('v', obj.v, 'residual', residual, 'u_t', u_t, 'f_x', f_x, 'viscosity', obj.RV.getViscosity(), 't', obj.t);
+            [residual, u_t, grad_f] = obj.RV.getResidual();
+            state = struct('v', obj.v, 'residual', residual, 'u_t', u_t, 'grad_f', grad_f, 'viscosity', obj.RV.getViscosity(), 't', obj.t);
         end
 
         function obj = step(obj)
-            F = @(v,t) obj.F(v,t,obj.RV.getViscosity());
-            v_p = obj.v;
-            obj.v = time.rk.rungekutta(obj.v, obj.t, obj.k, F, obj.coeffs);
+            obj.v = time.rk.rungekuttaRV(obj.v, obj.t, obj.k, obj.F, obj.RV, obj.coeffs);
             obj.t = obj.t + obj.k;
             obj.n = obj.n + 1;
-            obj.RV.update(obj.v,v_p,obj.k);
         end
     end