--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+time/+rk/General.m	Wed Jan 09 10:59:38 2019 +0100
@@ -0,0 +1,81 @@
+classdef General < time.Timestepper
+    properties
+        F       % RHS of the ODE
+        dt      % Time step
+        t       % Time point
+        v       % Solution vector
+        n       % Time level
+        scheme  % The scheme used for the time stepping, e.g rk4, rk6 etc.
+        coeffs  % Butcher tableau coefficients
+        V       % All stage approximations in most recent time step
+        K       % All stage rates in most recent time step
+    end
+
+
+    methods
+        % Timesteps v_t = F(v,t), using the specified RK method from t = t0 with
+        % timestep dt and initial conditions v = v0
+        function obj = General(F, dt, t0, v0, method, discreteData)
+            default_arg('method',"rk4");
+            default_arg('discreteData', []);
+            obj.F = F;
+            obj.dt = dt;
+            obj.t = t0;
+            obj.v = v0;
+            obj.n = 0;
+
+            % Extract the coefficients for the specified method
+            % used for the RK updates from the Butcher tableua.
+            [s,a,b,c] = time.rk.butcherTableauFromStr(method);
+            obj.coeffs = struct('s',s,'a',a,'b',b,'c',c);
+
+            if isempty(discreteData)
+                obj.scheme = @(v,t,n) time.rk.rungekutta(v, t, dt, F, obj.coeffs);
+            else
+                obj.scheme = @(v,t,n) time.rk.rungekuttaDiscreteData(v, t, dt, F, obj.coeffs, discreteData, n);
+            end
+        end
+
+        % v: Current solution
+        % t: Current time
+        % V: All stage approximations in most recent time step
+        % K: All stage rates in most recent time step
+        % T: Time points (corresponding to V and K) in most recent time step
+        function [v,t,V,T,K] = getV(obj)
+            v = obj.v;
+            t = obj.t;
+            V = obj.V;
+            K = obj.K;
+            T = obj.t + obj.dt*obj.coeffs.b;
+        end
+
+        function obj = step(obj)
+            [obj.v, obj.V, obj.K] = obj.scheme(obj.v, obj.t, obj.n);
+            obj.t = obj.t + obj.dt;
+            obj.n = obj.n + 1;
+        end
+
+        % Returns a vector of time points, including substage points,
+        % in the time interval [t0, tEnd].
+        % The time-step obj.dt is assumed to be aligned with [t0, tEnd] already.
+        function tvec = timePoints(obj, t0, tEnd)
+            N = round( (tEnd-t0)/obj.dt );
+            tvec = zeros(N*obj.s, 1);
+            s = obj.coeffs.s;
+            c = obj.coeffs.c;
+            for i = 1:N
+                ind = (i-1)*s+1 : i*s;
+                tvec(ind) = ((i-1) + c')*obj.dt;
+            end
+        end
+
+        % Returns a vector of quadrature weights corresponding to grid points
+        % in time interval [t0, tEnd], substage points included.
+        % The time-step obj.dt is assumed to be aligned with [t0, tEnd] already.
+        function weights = quadWeights(obj, t0, tEnd)
+            N = round( (tEnd-t0)/obj.dt );
+            b = obj.coeffs.b;
+            weights = repmat(b', N, 1);
+        end
+    end
+end
\ No newline at end of file