--- a/+time/Rungekutta4SecondOrder.m	Mon Sep 10 16:19:16 2018 +0200
+++ b/+time/Rungekutta4SecondOrder.m	Thu Nov 15 16:36:21 2018 -0800
@@ -15,6 +15,21 @@
 
 
     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;
@@ -22,21 +37,55 @@
             obj.m = length(v0);
             obj.n = 0;
 
-            I = speye(obj.m);
-            O = sparse(obj.m,obj.m);
-            obj.M = [O, I; D, E*I]; % Multiply with I to allow 0 as input.
+
+            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 S == 0
-                obj.C = zeros(2*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
-                obj.C = [zeros(obj.m,1), S];
-            end
+
+                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.k = k;
-            obj.t = t0;
-            obj.w = [v0; v0t];
+                obj.M = [
+                    O, I;
+                    D, E;
+                ];
+                obj.C = [
+                    zeros(obj.m,1);
+                                 S;
+                ];
 
-            obj.F = @(w,t)(obj.M*w + obj.C);
+                obj.k = k;
+                obj.t = t0;
+                obj.w = [v0; v0t];
+
+                obj.F = @(w,t)(obj.M*w + obj.C);
+            end
         end
 
         function [v,t] = getV(obj)