--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+sbp/+implementations/d1_gauss_4.m	Thu Feb 02 17:05:43 2017 +0100
@@ -0,0 +1,65 @@
+function [D1,H,x,h,e_l,e_r] = d1_gauss_4(N,L)
+
+% L: Domain length
+% N: Number of grid points
+if(nargin < 2)
+    L = 1;
+end
+
+if(N~=4)
+    error('This operator requires exactly 4 grid points');
+end
+
+% 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	Thu Feb 02 17:05:43 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(m,L);
+                otherwise
+                    error('Invalid operator order %d.',order);
+            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