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