--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+sbp/+implementations/d1_gauss_4.m	Tue Feb 07 15:59:22 2017 +0100
@@ -0,0 +1,63 @@
+function [D1,H,x,h,e_l,e_r] = d1_gauss_4(L)
+
+% L: Domain length
+% N: Number of grid points
+if(nargin < 2)
+    L = 1;
+end
+
+N = 4;
+
+% Quadrature nodes on interval [-1, 1]
+x = [ -0.8611363115940526; -0.3399810435848563; 0.3399810435848563; 0.8611363115940526];
+
+% Shift nodes to [0,L]
+x = (x+1)/2*L;
+
+% Boundary extrapolation operators
+e_l = [1.5267881254572668; -0.8136324494869273; 0.4007615203116504; -0.1139171962819899];
+e_r = flipud(e_l);
+e_l = sparse(e_l);
+e_r = sparse(e_r);
+
+%%%% Compute approximate h %%%%%%%%%%
+h = L/(N-1);
+%%%%%%%%%%%%%%%%%%%%%%%%%
+
+%%%% Norm matrix on [-1,1] %%%%%%%%
+P = sparse(N,N);
+P(1,1) =  0.3478548451374539;
+P(2,2) =  0.6521451548625461;
+P(3,3) =  0.6521451548625461;
+P(4,4) =  0.3478548451374539;
+%%%%%%%%%%%%%%%%%%%%%%%%%
+
+%%%% Norm matrix on [0,L] %%%%%%%%
+H = P*L/2;
+%%%%%%%%%%%%%%%%%%%%%%%%%
+
+%%%% D1 on [-1,1] %%%%%%%%
+D1 = sparse(N,N);
+D1(1,1) = -3.3320002363522817;
+D1(1,2) = 4.8601544156851962;
+D1(1,3) = -2.1087823484951789;
+D1(1,4) = 0.5806281691622644;
+
+D1(2,1) = -0.7575576147992339;
+D1(2,2) = -0.3844143922232086;
+D1(2,3) = 1.4706702312807167;
+D1(2,4) = -0.3286982242582743;
+
+D1(3,1) = 0.3286982242582743;
+D1(3,2) = -1.4706702312807167;
+D1(3,3) = 0.3844143922232086;
+D1(3,4) = 0.7575576147992339;
+
+D1(4,1) = -0.5806281691622644;
+D1(4,2) = 2.1087823484951789;
+D1(4,3) = -4.8601544156851962;
+D1(4,4) = 3.3320002363522817;
+%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% D1 on [0,L]
+D1 = D1*2/L;
\ No newline at end of file

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+sbp/D1Gauss.m	Tue Feb 07 15:59:22 2017 +0100
@@ -0,0 +1,41 @@
+classdef D1Gauss < sbp.OpSet
+    % Diagonal-norm SBP operators based on the Gauss quadrature formula
+    % with m nodes, which is of degree 2m-1. Hence, The operator D1 is
+    % accurate of order m.
+    properties
+        D1 % SBP operator approximating first derivative
+        H % Norm matrix
+        HI % H^-1
+        Q % Skew-symmetric matrix
+        e_l % Left boundary operator
+        e_r % Right boundary operator
+        m % Number of grid points.
+        h % Step size
+        x % grid
+        borrowing % Struct with borrowing limits for different norm matrices
+    end
+
+    methods
+        function obj = D1Gauss(m,lim)
+
+            x_l = lim{1};
+            x_r = lim{2};
+            L = x_r-x_l;
+
+            switch m
+                case 4
+                    [obj.D1,obj.H,obj.x,obj.h,obj.e_l,obj.e_r] = ...
+                        sbp.implementations.d1_gauss_4(L);
+                otherwise
+                    error('Invalid number of points: %d.', m);
+            end
+
+
+            obj.x = obj.x + x_l;
+            obj.HI = inv(obj.H);
+            obj.Q = obj.H*obj.D1 - obj.e_r*obj.e_r' + obj.e_l*obj.e_l';
+
+            obj.borrowing = [];
+        end
+    end
+end

--- a/+time/SBPInTime.m	Thu Feb 02 10:10:02 2017 +0100
+++ b/+time/SBPInTime.m	Tue Feb 07 15:59:22 2017 +0100
@@ -24,9 +24,16 @@
 
     methods
         function obj = SBPInTime(A, f, k, t0, v0, TYPE, order, blockSize)
-            default_arg('TYPE','minimal');
-            default_arg('order', 8);
-            default_arg('blockSize',time.SBPInTime.smallestBlockSize(order,TYPE));
+
+            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.SBPInTime.smallestBlockSize(order,TYPE));
+            end
 
             obj.A = A;
             obj.f = f;
@@ -45,6 +52,8 @@
                     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
 
             D1 = ops.D1;
@@ -76,7 +85,7 @@
 
             % Pretend that the initial condition is the last level
             % of a previous step.
-            obj.v = obj.Et_r * v0;
+            obj.v =  1/(e_r'*e_r) * obj.Et_r * v0;
 
         end
 
@@ -98,7 +107,7 @@
 
     methods(Static)
         function N = smallestBlockSize(order,TYPE)
-            default_arg('TYPE','equidistant')
+            default_arg('TYPE','gauss')
 
             switch TYPE
 
@@ -153,6 +162,8 @@
                         otherwise
                             error('Operator does not exist');
                     end
+                case 'gauss'
+                    N = 4;
             end
         end
     end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/Map.m	Tue Feb 07 15:59:22 2017 +0100
@@ -0,0 +1,91 @@
+classdef Map < handle
+    properties
+        map
+    end
+
+    % can we support multi map using varargin?
+    % probably a bad idea. For example it complicates keys();
+
+    methods
+        function obj = Map()
+            obj.map = containers.Map();
+        end
+
+        function set(obj, k, v)
+            keyByteStream = getByteStreamFromArray(k);
+
+            obj.map(char(keyByteStream)) = v;
+        end
+
+        function v = get(obj, k)
+            keyByteStream = getByteStreamFromArray(k);
+
+            v = obj.map(char(keyByteStream));
+        end
+
+        function b = isKey(obj, k)
+            keyByteStream = getByteStreamFromArray(k);
+            b = obj.map.isKey(char(keyByteStream));
+        end
+
+        function c = keys(obj)
+            keyByteStreams = obj.map.keys;
+
+            n = length(keyByteStreams);
+
+            c = cell(1, n);
+            for i = 1:n
+                c{i} = getArrayFromByteStream(uint8(keyByteStreams{i}));
+            end
+        end
+
+        function l = length(obj)
+            l = obj.map.length;
+        end
+
+        function remove(obj, k)
+            keyByteStream = getByteStreamFromArray(k);
+            obj.map.remove(char(keyByteStream));
+        end
+
+        function s = size(obj)
+            s = obj.map.size;
+        end
+
+        function c = values(obj)
+            c = obj.map.values;
+        end
+
+        function v = subsref(obj, S)
+            switch S(1).type
+                case '()'
+                    k = S.subs{1};
+                    try
+                        v = get(obj, k);
+                    catch ME
+                        if strcmp(ME.identifier,'MATLAB:Containers:Map:NoKey')
+                            error('Reference to non-existent entry %s',toString(S.subs));
+                        else
+                            throw(ME);
+                        end
+                    end
+                otherwise
+                    try
+                        v = builtin('subsref', obj, S);
+                    catch e
+                        error('You can''t use dot notation for this because Matlab(TM). What is this piece of shit software anyway?')
+                    end
+            end
+        end
+
+        function obj = subsasgn(obj, S, v);
+            switch S(1).type
+                case '()'
+                    k = S.subs{1};
+                    set(obj, k, v);
+                otherwise
+                    error('You can''t use dot notation because Matlab(TM). What is this piece of shit software anyway?')
+            end
+        end
+    end
+end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/MapTest.m	Tue Feb 07 15:59:22 2017 +0100
@@ -0,0 +1,97 @@
+function tests = MatTest()
+    tests = functiontests(localfunctions);
+end
+
+function kvp = getKeyValuePairs()
+    kvp = {
+        {1,3},1;
+        struct(), [1; 3; 4];
+        [1,2; 4 3], struct();
+        'Hej', struct('lol', 6);
+        0, 'Nej';
+    };
+end
+
+function testSetAndGet(testCase)
+    keyValuePairs = getKeyValuePairs();
+
+    map = Map();
+
+    % Insert
+    for i = 1:length(keyValuePairs)
+        map(keyValuePairs{i,1}) = keyValuePairs{i,2};
+    end
+
+    % Validate output
+    for i = 1:length(keyValuePairs)
+        v = map(keyValuePairs{i,1});
+        testCase.verifyEqual(v, keyValuePairs{i,2});
+    end
+end
+
+function map = exampleMap()
+    keyValuePairs = getKeyValuePairs();
+
+    map = Map();
+
+    % Insert
+    for i = 1:length(keyValuePairs)
+        map(keyValuePairs{i,1}) = keyValuePairs{i,2};
+    end
+end
+
+function testLength(testCase)
+    map = Map();
+    testCase.verifyEqual(map.length, 0);
+
+    map = exampleMap();
+    testCase.verifyEqual(map.length, 5)
+end
+
+
+function testIsKey(testCase)
+    map = exampleMap();
+
+    keyValuePairs = getKeyValuePairs();
+    keys = keyValuePairs(:,1);
+
+    for i = 1:length(keys)
+        testCase.verifyTrue(map.isKey(keys{i}));
+    end
+
+    notKeys = {
+        'hej',
+        [],
+        1,
+        {2,5},
+    };
+
+    for i = 1:length(notKeys)
+        testCase.verifyFalse(map.isKey(notKeys{i}));
+    end
+end
+
+
+function testRemove(testCase)
+    map = exampleMap();
+
+    remove(map, struct());
+
+    testCase.verifyFalse(map.isKey(struct()));
+end
+
+% function testValues(testCase)
+%     keyValuePairs = getKeyValuePairs();
+
+%     map = exampleMap();
+
+%     testCase.verifyEqual(values(map), keyValuePairs(:,2)');
+% end
+
+% function testKeys(testCase)
+%     keyValuePairs = getKeyValuePairs();
+
+%     map = exampleMap();
+
+%     testCase.verifyEqual(keys(map), keyValuePairs(:,1)');
+% end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/diags.m	Tue Feb 07 15:59:22 2017 +0100
@@ -0,0 +1,9 @@
+function A = diags(B,d,m,n)
+    assert(size(B,1) == m);
+
+    A = repmat(B(:,1)*0, [1, n]);
+
+    for i = 1:size(B,2)
+        A(:,d(i)+ (1:m)) = A(:,d(i)+ (1:m)) + diag(B(:,i));
+    end
+end