--- a/diracDiscr.m	Wed Nov 15 14:56:52 2017 -0800
+++ b/diracDiscr.m	Wed Nov 15 14:58:13 2017 -0800
@@ -1,5 +1,7 @@
 function ret = diracDiscr(x_0in , x , m_order, s_order, H)
 
+m = length(x);
+
 % Return zeros if x0 is outside grid
 if(x_0in < x(1) || x_0in > x(end) )
 
@@ -14,7 +16,7 @@
     S = [];
     M = [];
     poss = find(tot*h/2 >= eta);
-
+    
     % Ensure that poss is not too long
     if length(poss) == (tot + 2)
         poss = poss(2:end-1);
@@ -23,19 +25,31 @@
     end
 
     % Use first tot grid points
-    if length(poss)<tot && eta(end)>eta(1)
+    if length(poss)<tot && x_0in < x(1) + ceil(tot/2)*h;
         index=1:tot;
         pol=(x(1:tot)-x(1))/(x(tot)-x(1));
         x_0=(x_0in-x(1))/(x(tot)-x(1));
         norm=fnorm(1:tot)/h;
 
     % Use last tot grid points
-    elseif length(poss)<tot && eta(end)<eta(1)
+    elseif length(poss)<tot && x_0in > x(end) - ceil(tot/2)*h;
         index = length(x)-tot+1:length(x);
         pol = (x(end-tot+1:end)-x(end-tot+1))/(x(end)-x(end-tot+1));
         norm = fnorm(end-tot+1:end)/h;
         x_0 = (x_0in-x(end-tot+1))/(x(end)-x(end-tot+1));
 
+    % Interior, compensate for round-off errors.
+    elseif length(poss) < tot
+        if poss(end)<m
+            poss = [poss; poss(end)+1];
+        else
+            poss = [poss(1)-1; poss];
+        end
+        pol = (x(poss)-x(poss(1)))/(x(poss(end))-x(poss(1)));
+        x_0 = (x_0in-x(poss(1)))/(x(poss(end))-x(poss(1)));
+        norm = fnorm(poss)/h;
+        index = poss;
+        
     % Interior
     else    
         pol = (x(poss)-x(poss(1)))/(x(poss(end))-x(poss(1)));

--- a/diracDiscrTest.m	Wed Nov 15 14:56:52 2017 -0800
+++ b/diracDiscrTest.m	Wed Nov 15 14:58:13 2017 -0800
@@ -193,12 +193,39 @@
     end
 end
 
+function testAllGP(testCase)
+
+    orders = [2, 4, 6];
+    mom_conds = orders;
+    
+    for o = 1:length(orders)
+        order = orders(o);
+        mom_cond = mom_conds(o);
+        [xl, xr, m, h, x, H, fs] = setupStuff(order, mom_cond);
+
+        % Test all grid points
+        x0s = x;
+
+        for j = 1:length(fs)
+                f = fs{j};
+                fx = f(x);
+            for i = 1:length(x0s)
+                x0 = x0s(i);
+                delta = diracDiscr(x0, x, mom_cond, 0, H);
+                integral = delta'*H*fx;
+                err = abs(integral - f(x0));
+                testCase.verifyLessThan(err, 1e-12);
+            end
+        end
+    end
+end
+
 function [xl, xr, m, h, x, H, fs] = setupStuff(order, mom_cond)
 
     % Grid
     xl = -3;
     xr = 2;
-    m = 21;
+    m = 101;
     h = (xr-xl)/(m-1);
     g = grid.equidistant(m, {xl, xr});
     x = g.points();