--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+rv/+time/BdfDerivative.m	Fri Dec 07 13:11:53 2018 +0100
@@ -0,0 +1,31 @@
+classdef BdfDerivative < handle
+    properties
+        coefficients
+    end
+    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
+        % 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