diff -r 1147db8a2ffa -r e0caae9ef6ed +time/SBPInTimeImplicitFormulation.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+time/SBPInTimeImplicitFormulation.m	Mon Jun 19 16:50:13 2017 +0200
@@ -0,0 +1,124 @@
+classdef SBPInTimeImplicitFormulation < time.Timestepper
+    % The SBP in time method.
+    % Implemented for A*v_t = B*v + f(t), v(0) = v0
+    properties
+        A,B
+        f
+
+        k % total time step.
+
+        blockSize % number of points in each block
+        N % Number of components
+
+        order
+        nodes
+
+        M,K     % System matrices
+        L,U,p,q % LU factorization of M
+        e_T
+
+        % Time state
+        t
+        v
+        n
+    end
+
+    methods
+        function obj = SBPInTimeImplicitFormulation(A, B, f, k, t0, v0, TYPE, order, blockSize)
+
+            default_arg('TYPE','gauss');
+
+            if(strcmp(TYPE,'gauss'))
+                default_arg('order',4)
+                default_arg('blockSize',4)
+            else
+                default_arg('order', 8);
+                default_arg('blockSize',time.SBPInTimeImplicitFormulation.smallestBlockSize(order,TYPE));
+            end
+
+            obj.A = A;
+            obj.B = B;
+            obj.f = f;
+
+            obj.k = k;
+            obj.blockSize = blockSize;
+            obj.N = length(v0);
+
+            obj.n = 0;
+            obj.t = t0;
+
+            %==== Build the time discretization matrix =====%
+            switch TYPE
+                case 'equidistant'
+                    ops = sbp.D2Standard(blockSize,{0,obj.k},order);
+                case 'optimal'
+                    ops = sbp.D1Nonequidistant(blockSize,{0,obj.k},order);
+                case 'minimal'
+                    ops = sbp.D1Nonequidistant(blockSize,{0,obj.k},order,'minimal');
+                case 'gauss'
+                    ops = sbp.D1Gauss(blockSize,{0,obj.k});
+            end
+
+            I = speye(size(A));
+            I_t = speye(blockSize,blockSize);
+
+            D1 = kron(ops.D1, I);
+            HI = kron(ops.HI, I);
+            e_0 = kron(ops.e_l, I);
+            e_T = kron(ops.e_r, I);
+            obj.nodes = ops.x;
+
+            % Convert to form M*w = K*v0 + f(t)
+            tau = kron(I_t, A) * e_0;
+            M = kron(I_t, A)*D1 + HI*tau*e_0' - kron(I_t, B);
+
+            K = HI*tau;
+
+            obj.M = M;
+            obj.K = K;
+            obj.e_T = e_T;
+
+            % LU factorization
+            [obj.L,obj.U,obj.p,obj.q] = lu(obj.M, 'vector');
+
+            obj.v = v0;
+        end
+
+        function [v,t] = getV(obj)
+            v = obj.v;
+            t = obj.t;
+        end
+
+        function obj = step(obj)
+            RHS = zeros(obj.blockSize*obj.N,1);
+
+            for i = 1:length(obj.blockSize)
+                RHS((1 + (i-1)*obj.N):(i*obj.N)) = obj.f(obj.nodes(i));
+            end
+
+            RHS = RHS + obj.K*obj.v;
+
+            y = obj.L\RHS(obj.p);
+            z = obj.U\y;
+
+            w = zeros(size(z));
+            w(obj.q) = z;
+
+            obj.v = obj.e_T'*w;
+
+            obj.t = obj.t + obj.k;
+            obj.n = obj.n + 1;
+        end
+    end
+
+    methods(Static)
+        function N = smallestBlockSize(order,TYPE)
+            default_arg('TYPE','gauss')
+
+            switch TYPE
+                case 'gauss'
+                    N = 4;
+            end
+        end
+    end
+end

diff -r 1147db8a2ffa -r e0caae9ef6ed +time/SBPInTimeSecondOrderFormImplicit.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+time/SBPInTimeSecondOrderFormImplicit.m	Mon Jun 19 16:50:13 2017 +0200
@@ -0,0 +1,80 @@
+classdef SBPInTimeSecondOrderFormImplicit < time.Timestepper
+    properties
+        A, B, C, f
+        AA, BB, ff
+
+        n
+        t
+        k
+
+        firstOrderTimeStepper
+    end
+
+    methods
+        % Solves A*u_tt + B*u_t + C*u + f(t) = 0
+        % A, B can either both be constants or both be function handles,
+        % They can also be omitted by setting them equal to the empty matrix.
+        function obj = SBPInTimeSecondOrderFormImplicit(A, B, C, f, k, t0, v0, v0t, TYPE, order, blockSize)
+            default_arg('f', []);
+            default_arg('TYPE', []);
+            default_arg('order', []);
+            default_arg('blockSize',[]);
+
+            m = length(v0);
+
+            default_arg('A', sparse(m, m));
+            default_arg('B', sparse(m, m));
+            default_arg('C', sparse(m, m));
+
+            I = speye(m);
+            O = sparse(m,m);
+
+            % Rewrite to
+            % AA*w_t = BB*w + ff(t);
+
+            obj.AA = [
+                 I, O;
+                 O, A;
+            ];
+            obj.BB = [
+                 O,  I;
+                -B, -C;
+            ];
+
+            if ~isempty(f)
+                obj.ff = @(t)[
+                    sparse(m,1);
+                           -f(t);
+                ];
+            else
+                obj.ff = @(t) sparse(2*m,1);
+            end
+
+            w0 = [v0; v0t];
+
+            obj.k = k;
+            obj.t = t0;
+            obj.n = 0;
+
+            obj.firstOrderTimeStepper = time.SBPInTimeImplicitFormulation(obj.AA, obj.BB, obj.ff, obj.k, obj.t, w0, TYPE, order, blockSize);
+        end
+
+        function [v,t] = getV(obj)
+            w = obj.firstOrderTimeStepper.getV();
+            v = w(1:end/2);
+            t = obj.t;
+        end
+
+        function [vt,t] = getVt(obj)
+            w = obj.firstOrderTimeStepper.getV();
+            vt = w(end/2+1:end);
+            t = obj.t;
+        end
+
+        function obj = step(obj)
+            obj.firstOrderTimeStepper.step();
+            obj.t = obj.firstOrderTimeStepper.t;
+            obj.n = obj.firstOrderTimeStepper.n;
+        end
+    end
+end