--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+sbp/+implementations/d1_gauss_4.m	Tue Feb 07 12:54:52 2017 +0000
@@ -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 12:54:52 2017 +0000
@@ -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 14:46:36 2017 +0000
+++ b/+time/SBPInTime.m	Tue Feb 07 12:54:52 2017 +0000
@@ -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
 
@@ -97,7 +106,7 @@
 
     methods(Static)
         function N = smallestBlockSize(order,TYPE)
-            default_arg('TYPE','equidistant')
+            default_arg('TYPE','gauss')
 
             switch TYPE
 
@@ -152,6 +161,8 @@
                         otherwise
                             error('Operator does not exist');
                     end
+                case 'gauss'
+                    N = 4;
             end
         end
     end