--- a/+rv/+time/BDFDerivative.m	Wed Dec 05 15:04:44 2018 +0100
+++ b/+rv/+time/BDFDerivative.m	Thu Dec 06 11:30:47 2018 +0100
@@ -1,33 +1,31 @@
-class BDFDerivative < handle
+classdef BDFDerivative < handle
     properties
-        dt
-        coefficents
+        coefficients
     end
-    methods
-        function obj = BDFDerivative(order, dt)
-            obj.coefficents = obj.getBDFCoefficients(order);
-            obj.dt = dt;
-        end
-
-        function DvDt = evaluate(v, v_prev)
-            DvDt = (obj.coefficients(1)*v - sum(obj.coefficients(2:end).*v_prev),2)/obj.dt;
+    methods %TBD: Decide if the BDF order should be passed in construction
+            %     and only a row of coefficients stored based on the order.
+            %     There would still be an implicit dependancy on the number
+            %     of vectors in v_prev and elements in coefficients.
+            %     In addition, dt could be stored, but this would be
+            %     inflexible if different step sizes are employed.
+        function obj = BDFDerivative()
+            obj.coefficients = obj.getBDFCoefficients();
         end
-
-        function c = getBDFCoefficients(obj, order)
-            switch order
-                case 1
-                    c(1,1) = 1; c(1,2) = 1;
-                case 2 
-                    c(1,1) = 3/2; c(1,2) = 4/2; c(1,3) = -1/2;
-                case 3
-                    c(1,1) = 11/6; c(1,2) = 18/6; c(1,3) = -9/6; c(1,4) = 2/6;
-                case 4
-                    c(1,1) = 25/12; c(1,2) = 48/12; c(1,3) = -36/12; c(1,4) = 16/12; c(1,5) = -3/12;
-                case 5
-                    c(1,1) = 137/60; c(1,2) = 300/60; c(1,3) = -300/60; c(1,4) = 200/60; c(1,5) = -75/60; c(1,6) = 12/60;
-                case 6
-                    c(1,1) = 147/60; c(1,2) = 360/60; c(1,3) = -450/60; c(1,4) = 400/60; c(1,5) = -225/60; c(1,6) = 72/60; c(1,7) = -10/60;
-            assert(length(c) == (order+1));
+        % Add asserts here?
+        function DvDt = evaluate(obj, v, v_prev, dt)
+            order = size(v_prev,2);
+            DvDt = (obj.coefficients(order,1)*v - sum(obj.coefficients(order,2:order+1).*v_prev,2))/dt;
+        end
+    end
+    methods(Static)
+        function c = getBDFCoefficients()
+            c = zeros(6,7);
+            c(1,1) = 1; c(1,2) = 1;
+            c(2,1) = 3/2; c(2,2) = 4/2; c(2,3) = -1/2;
+            c(3,1) = 11/6; c(3,2) = 18/6; c(3,3) = -9/6; c(3,4) = 2/6;
+            c(4,1) = 25/12; c(4,2) = 48/12; c(4,3) = -36/12; c(4,4) = 16/12; c(4,5) = -3/12;
+            c(5,1) = 137/60; c(5,2) = 300/60; c(5,3) = -300/60; c(5,4) = 200/60; c(5,5) = -75/60; c(5,6) = 12/60;
+            c(6,1) = 147/60; c(6,2) = 360/60; c(6,3) = -450/60; c(6,4) = 400/60; c(6,5) = -225/60; c(6,6) = 72/60; c(6,7) = -10/60;
         end
     end
 end

--- a/+rv/+time/RungekuttaExteriorRV.m	Wed Dec 05 15:04:44 2018 +0100
+++ b/+rv/+time/RungekuttaExteriorRV.m	Thu Dec 06 11:30:47 2018 +0100
@@ -7,6 +7,7 @@
         n       % Time level
         coeffs  % The coefficents used for the RK time integration
         RV      % Residual Viscosity
+        bdfOrder
         v_prev  % Solution vector at previous time levels, used for the RV update
         DvDt    % Function for computing the time deriative used for the RV update
     end
@@ -24,10 +25,14 @@
             obj.coeffs = struct('s',s,'a',a,'b',b,'c',c);
         
             obj.RV = RV;
-            % TODO: Initialize v_prev.
-            obj.v_prev = zeros(length(obj.v),bdfOrder);
             %  TBD: For cases where h~k we could probably use rkOrder-2 here.
-            obj.DvDt = rv.time.BDFDerivative(bdfOrder,k);
+            %  TBD: Decide on if the initialization of the previous stages used by
+            %       the BDF should be done here, or if it should be checked for each
+            %       step taken.
+            assert(bdfOrder >= 2);
+            obj.bdfOrder = bdfOrder;
+            obj.v_prev = [];
+            obj.DvDt = rv.time.BDFDerivative();
         end
 
         function [v, t] = getV(obj)
@@ -41,15 +46,19 @@
         end
 
         function obj = step(obj)
-            % Store current time level    
-            obj.v_prev(:,2:end) = obj.v_prev(:,1:end-1);
-            obj.v_prev(:,1) = obj.v;       
+            % Store current time level
+            if (size(obj.v_prev,2) < obj.bdfOrder)
+                obj.v_prev = [obj.v, obj.v_prev];
+            else
+                obj.v_prev(:,2:end) = obj.v_prev(:,1:end-1);
+                obj.v_prev(:,1) = obj.v;
+            end    
             % Fix the viscosity of the RHS function F
-            F_visc = @(v,t)F(v,t,obj.RV.getViscosity()); %Remove state in RV?
-            obj.v = time.rk.rungekutta(obj.v, obj.t, obj.k, F_visc, obj.RV, obj.coeffs);
+            F_visc = @(v,t) obj.F(v,t,obj.RV.getViscosity()); %TBD: Remove state in RV?
+            obj.v = time.rk.rungekutta(obj.v, obj.t, obj.k, F_visc, obj.coeffs);
             % Calculate dvdt and update RV
-            dvdt = obj.DvDt.evaluate(v, v_prev, lenght(v_prev));
-            RV.update(v,dvdt);
+            dvdt = obj.DvDt.evaluate(obj.v, obj.v_prev, obj.k);
+            obj.RV.update(obj.v,dvdt);
             obj.t = obj.t + obj.k;
             obj.n = obj.n + 1;
         end

--- a/+rv/ResidualViscosity.m	Wed Dec 05 15:04:44 2018 +0100
+++ b/+rv/ResidualViscosity.m	Thu Dec 06 11:30:47 2018 +0100
@@ -40,6 +40,13 @@
             obj.viscosity = min(obj.Cmax*obj.h*abs(obj.waveSpeed(v)), obj.Cres*obj.h^2*abs(obj.residual)/norm(v-mean(v),inf));
         end
 
+        function smoothendVector = smoothen(obj, vector)
+            smoothendVector = vector;
+            for i = 2:length(vector)-1
+                smoothendVector(i) = (1/6)*(vector(i-1) + 4*vector(i) + vector(i+1));
+            end
+        end
+
         function [residual, u_t, grad_f] = getResidual(obj)
             residual = obj.residual;
             u_t = obj.u_t;