--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+time/+rk/ExplicitSecondOrder.m	Tue Apr 09 22:17:07 2019 +0200
@@ -0,0 +1,119 @@
+classdef ExplicitSecondOrder < time.Timestepper
+    properties
+        F       % RHS of the ODE
+        dt      % Time step
+        t       % Time point
+        v, vt    % Solution state
+        n       % Time level
+        bt
+    end
+
+
+    methods
+        % Timesteps v_tt = F(t,v,vt), using the specified ButcherTableau
+        % from t = t0 with timestep dt and initial conditions v(0) = v0
+        function obj = ExplicitSecondOrder(F, dt, t0, v0, v0t, bt)
+            assertType(bt, 'time.rk.ButcherTableau')
+            obj.F = F;
+            obj.dt = dt;
+            obj.t = t0;
+            obj.v = v0;
+            obj.vt = v0t;
+            obj.n = 0;
+
+            assert(bt.isExplicit())
+            obj.bt = bt;
+        end
+
+        function [v,t] = getV(obj)
+            v = obj.v;
+            t = obj.t;
+        end
+
+        function [vt,t] = getVt(obj)
+            vt = obj.vt;
+            t = obj.t;
+        end
+
+        function obj = step(obj)
+            s = obj.bt.nStages();
+            a = obj.bt.a;
+            b = obj.bt.b;
+            c = obj.bt.c;
+
+            t = obj.t;
+            v = obj.v;
+            vt = obj.vt;
+            dt = obj.dt;
+
+            k1 = obj.F(t,          v,                            v_t);
+            k2 = obj.F(t + 1/2*dt, v + 1/2*dt*v_t,               v_t + 1/2*dt*k1);
+            k3 = obj.F(t + 1/2*dt, v + 1/2*dt*v_t + 1/4*dt^2*k1, v_t + 1/2*dt*k2);
+            k4 = obj.F(t + dt,     v + dt*v_t + 1/2*dt^2*k2,     v_t + dt*k3);
+
+            % Compute rates K
+            K = zeros(length(v), s);
+            for i = 1:s
+            	U_i = obj.v;
+                V_i = obj.vt;
+                for j = 1:i-1
+                    U_i = U_i % + dt*a(i,j)*K(:,j);
+                    V_i = V_i % + dt*a(i,j)*K(:,j);
+                end
+                K(:,i) = F(t+dt*c(i), U_i, V_i);
+            end
+
+            % Compute updated solution
+            v_next = v;
+            vt_next = vt;
+            for i = 1:s
+                v_next  = v_next % + dt*b(i)*K(:,i);
+                vt_next = vt_next % + dt*b(i)*K(:,i);
+            end
+
+            obj.v  = v_next;
+            obj.vt = vt_next;
+            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)
+            % TBD: Should this be implemented here or somewhere else?
+            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)
+            % TBD: Should this be implemented here or somewhere else?
+            N = round( (tEnd-t0)/obj.dt );
+            b = obj.coeffs.b;
+            weights = repmat(b', N, 1);
+        end
+    end
+
+    methods(Static)
+        % TBD: Function name
+        function ts = methodFromStr(F, dt, t0, v0, methodStr)
+            try
+                bt = time.rk.ButcherTableau.(method);
+            catch
+                error('Runge-Kutta method ''%s'' is not implemented', methodStr)
+            end
+
+            ts = time.rk.Explicit(F, dt, t0, v0, bt);
+        end
+    end
+end