--- a/+time/Rk4SecondOrderNonlin.m	Tue Jan 08 12:19:33 2019 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,85 +0,0 @@
-classdef Rk4SecondOrderNonlin < time.Timestepper
-    properties
-        F
-        k
-        t
-        w
-        m
-
-        D
-        E
-        S
-
-        n
-    end
-
-
-    methods
-        function obj = Rk4SecondOrderNonlin(D, E, S, k, t0, v0, v0t)
-            default_arg('S',0);
-            default_arg('E',0);
-
-            if isnumeric(S)
-                S = @(v,t)S;
-            end
-
-            if isnumeric(E)
-                E = @(v)E;
-            end
-
-            obj.k = k;
-            obj.t = t0;
-            obj.w = [v0; v0t];
-
-            m = length(v0);
-            function wt = F(w,t)
-                v  = w(1:m);
-                vt = w(m+1:end);
-
-                % Def: w = [v; vt]
-                wt(1:m,1) = vt;
-                wt(m+1:2*m,1) = D(v)*v + E(v)*vt + S(v,t);
-
-            end
-
-            obj.F = @F;
-            obj.D = D;
-            obj.E = E;
-            obj.S = S;
-            obj.m = m;
-            obj.n = 0;
-        end
-
-        function [v,t] = getV(obj)
-            v = obj.w(1:end/2);
-            t = obj.t;
-        end
-
-        function [vt,t] = getVt(obj)
-            vt = obj.w(end/2+1:end);
-            t = obj.t;
-        end
-
-        function obj = step(obj)
-            w = obj.w;
-            k = obj.k;
-
-            k1 = obj.F(t, w);
-            k2 = obj.F(t + 0.5*k, w + 0.5*k*k1);
-            k3 = obj.F(t + 0.5*k, w + 0.5*k*k2);
-            k4 = obj.F(t + k, w + k*k3);
-
-            obj.w = w + k*(1/6)*(k1+2*(k2+k3)+k4);
-            obj.t = obj.t + obj.k;
-            obj.n = obj.n + 1;
-        end
-    end
-
-
-    methods (Static)
-        function k = getTimeStep(lambda)
-            k = rk4.get_rk4_time_step(lambda);
-        end
-    end
-
-end
\ No newline at end of file

--- a/+time/Rungekutta4SecondOrder.m	Tue Jan 08 12:19:33 2019 +0100
+++ b/+time/Rungekutta4SecondOrder.m	Tue Jan 08 12:34:29 2019 +0100
@@ -1,91 +1,31 @@
 classdef Rungekutta4SecondOrder < time.Timestepper
     properties
         F
-        k
+        dt
         t
+        n
+
+        G % RHS for rewritten system
         w
-        m
-        D
-        E
-        S
-        M
-        C
-        n
     end
 
 
     methods
-        % Solves u_tt = Du + Eu_t + S by
-        % Rewriting on first order form:
-        %   w_t = M*w + C(t)
-        % where
-        %   M = [
-        %      0, I;
-        %      D, E;
-        %   ]
-        % and
-        %   C(t) = [
-        %      0;
-        %      S(t)
-        %   ]
-        % D, E, S can either all be constants or all be function handles,
-        % They can also be omitted by setting them equal to the empty matrix.
-        function obj = Rungekutta4SecondOrder(D, E, S, k, t0, v0, v0t)
-            obj.D = D;
-            obj.E = E;
-            obj.S = S;
-            obj.m = length(v0);
+        % Create a time stepper for
+        %   v_tt = F(t,v,v_t),  v(t0) = v0, v_t(t0) = v0t
+        % with step size dt, by rewriting on first order form
+        function obj = Rungekutta4SecondOrder(F, dt, t0, v0, v0t)
+            obj.F = F;
+            obj.dt = dt;
+            obj.t = t0;
             obj.n = 0;
 
-
-            if isa(D, 'function_handle') || isa(E, 'function_handle') || isa(S, 'function_handle')
-                default_arg('D', @(t)sparse(obj.m, obj.m));
-                default_arg('E', @(t)sparse(obj.m, obj.m));
-                default_arg('S', @(t)sparse(obj.m, 1)    );
-
-                if ~isa(D, 'function_handle')
-                    D = @(t)D;
-                end
-                if ~isa(E, 'function_handle')
-                    E = @(t)E;
-                end
-                if ~isa(S, 'function_handle')
-                    S = @(t)S;
-                end
-
-                obj.k = k;
-                obj.t = t0;
-                obj.w = [v0; v0t];
-
-                % Avoid matrix formulation because it is VERY slow
-                obj.F = @(w,t)[
-                    w(obj.m+1:end);
-                    D(t)*w(1:obj.m) + E(t)*w(obj.m+1:end) + S(t);
-                ];
-            else
-
-                default_arg('D', sparse(obj.m, obj.m));
-                default_arg('E', sparse(obj.m, obj.m));
-                default_arg('S', sparse(obj.m, 1)    );
-
-                I = speye(obj.m);
-                O = sparse(obj.m,obj.m);
-
-                obj.M = [
-                    O, I;
-                    D, E;
-                ];
-                obj.C = [
-                    zeros(obj.m,1);
-                                 S;
-                ];
-
-                obj.k = k;
-                obj.t = t0;
-                obj.w = [v0; v0t];
-
-                obj.F = @(w,t)(obj.M*w + obj.C);
-            end
+            m = length(v0);
+            obj.w = [v0; v0t];
+            obj.G = @(t, w)[
+                w(m+1:end);
+                F(t,w(1:m), w(m+1:end));
+            ];
         end
 
         function [v,t] = getV(obj)
@@ -99,16 +39,17 @@
         end
 
         function obj = step(obj)
+            % TODO: Avoid extra storage
             w = obj.w;
-            k = obj.k;
+            dt = obj.dt;
 
-            k1 = obj.F(t, w);
-            k2 = obj.F(t + 0.5*k, w + 0.5*k*k1);
-            k3 = obj.F(t + 0.5*k, w + 0.5*k*k2);
-            k4 = obj.F(t + k, w + k*k3);
+            k1 = obj.G(t, w);
+            k2 = obj.G(t + 0.5*dt, w + 0.5*dt*k1);
+            k3 = obj.G(t + 0.5*dt, w + 0.5*dt*k2);
+            k4 = obj.G(t + dt, w + dt*k3);
 
-            obj.w = w + k*(1/6)*(k1+2*(k2+k3)+k4);
-            obj.t = obj.t + obj.k;
+            obj.w = w + dt*(1/6)*(k1+2*(k2+k3)+k4);
+            obj.t = obj.t + obj.dt;
             obj.n = obj.n + 1;
         end
     end