--- a/+time/SBPInTime.m	Wed Feb 08 09:07:07 2017 +0100
+++ b/+time/SBPInTime.m	Wed Feb 08 09:18:08 2017 +0100
@@ -1,89 +1,92 @@
 classdef SBPInTime < time.Timestepper
     % The SBP in time method.
     % Implemented for v_t = A*v + f(t)
-    % k_local -- time-step
-    % Nblock -- number of points in each block
-    % nodes -- points such that t_n + nodes are the points in block n.
-    % Each "step" takes one block step and thus advances 
-    % k = k_local*(Nblock-1) in time.
-    % M -- matrix used in every solve.
-    % [L,U,P,Q] = lu(M);
+    %
+    % Each "step" takes one block step and thus advances
+    % k = k_local*(blockSize-1) in time.
     properties
-        M
-        L
-        U
-        P
-        Q
+        M     % System matrix
+        L,U,P,Q % LU factorization of M
         A
         Et_r
         penalty
         f
-        k_local
-        k
+        k_local % step size within a block
+        k % Time size of a block  k/(blockSize-1) = k_local
         t
         v
         m
         n
-        Nblock
+        blockSize % number of points in each block
         order
         nodes
     end
 
     methods
-        function obj = SBPInTime(A, f, k, order, Nblock, t0, v0, TYPE)
-            default_arg('TYPE','equidistant');
-            default_arg('Nblock',time.SBPInTime.smallestBlockSize(order,TYPE));
-            
+        function obj = SBPInTime(A, 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.SBPInTime.smallestBlockSize(order,TYPE));
+            end
+
             obj.A = A;
             obj.f = f;
-            obj.k_local = k;
-            obj.k = k*(Nblock-1);
-            obj.Nblock = Nblock;
+            obj.k_local = k/(blockSize-1);
+            obj.k = k;
+            obj.blockSize = blockSize;
             obj.t = t0;
             obj.m = length(v0);
             obj.n = 0;
-            
+
             %==== Build the time discretization matrix =====%
             switch TYPE
                 case 'equidistant'
-                    ops = sbp.D2Standard(Nblock,{0,obj.k},order);
+                    ops = sbp.D2Standard(blockSize,{0,obj.k},order);
                 case 'optimal'
-                    ops = sbp.D1Nonequidistant(Nblock,{0,obj.k},order);
+                    ops = sbp.D1Nonequidistant(blockSize,{0,obj.k},order);
                 case 'minimal'
-                    ops = sbp.D1Nonequidistant(Nblock,{0,obj.k},order,'minimal');
+                    ops = sbp.D1Nonequidistant(blockSize,{0,obj.k},order,'minimal');
+                case 'gauss'
+                    ops = sbp.D1Gauss(blockSize,{0,obj.k});
             end
-            
+
             D1 = ops.D1;
             HI = ops.HI;
             e_l = ops.e_l;
             e_r = ops.e_r;
             obj.nodes = ops.x;
-            
+
             Ix = speye(size(A));
-            It = speye(Nblock,Nblock);
-            
-            obj.Et_r = kron(e_r,Ix);            
-            
+            It = speye(blockSize,blockSize);
+
+            obj.Et_r = kron(e_r,Ix);
+
             % Time derivative + penalty
             tau = 1;
-            Mt = D1 + tau*HI*(e_l*e_l'); 
-            
+            Mt = D1 + tau*HI*(e_l*e_l');
+
             % penalty to impose "data"
             penalty = tau*HI*e_l;
             obj.penalty = kron(penalty,Ix);
-            
+
             Mx = kron(It,A);
-            Mt = kron(Mt,Ix); 
+            Mt = kron(Mt,Ix);
             obj.M = Mt - Mx;
             %==============================================%
-            
+
             % LU factorization
             [obj.L,obj.U,obj.P,obj.Q] = lu(obj.M);
-            
+
             % 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
 
         function [v,t] = getV(obj)
@@ -93,39 +96,20 @@
 
         function obj = step(obj)
             obj.v = time.sbp.sbpintime(obj.v, obj.t, obj.nodes,...
-                              obj.penalty, obj.f, obj.Nblock,...
+                              obj.penalty, obj.f, obj.blockSize,...
                               obj.Et_r,...
                               obj.L, obj.U, obj.P, obj.Q);
             obj.t = obj.t + obj.k;
-            obj.n = obj.n + obj.Nblock-1;
+            obj.n = obj.n + 1;
         end
     end
-    
-    
+
     methods(Static)
-        
-        % 
-        function [k,numberOfBlocks] = alignedTimeStep(k,Tend,Nblock)
-            
-            % input k is the desired time-step
-            % Nblock is the number of points per block.
-            
-            % Make sure that we reach the final time by advancing 
-            % an integer number of blocks
-            kblock = (Nblock-1)*k;
-            numberOfBlocks = ceil(Tend/kblock);
-            kblock = Tend/(numberOfBlocks);
+        function N = smallestBlockSize(order,TYPE)
+            default_arg('TYPE','gauss')
 
-            % Corrected time step
-            k = kblock/(Nblock-1);
-            
-        end
-        
-        function N = smallestBlockSize(order,TYPE)
-            default_arg('TYPE','equidistant')
-            
             switch TYPE
-                
+
                 case 'equidistant'
                     switch order
                         case 2
@@ -143,9 +127,9 @@
                         otherwise
                             error('Operator does not exist');
                     end
-                    
+
                 case 'optimal'
-                    
+
                     switch order
                         case 4
                             N = 8;
@@ -160,9 +144,9 @@
                         otherwise
                             error('Operator does not exist');
                     end
-                    
+
                 case 'minimal'
-                    
+
                     switch order
                         case 4
                             N = 6;
@@ -177,13 +161,9 @@
                         otherwise
                             error('Operator does not exist');
                     end
-                    
+                case 'gauss'
+                    N = 4;
             end
-            
         end
-        
     end
-    
-    
-
 end