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