--- 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();