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