--- a/diracDiscr.m	Tue Nov 14 15:35:51 2017 -0800
+++ b/diracDiscr.m	Tue Nov 14 15:40:06 2017 -0800
@@ -1,66 +1,75 @@
 function ret = diracDiscr(x_0in , x , m_order, s_order, H)
 
-fnorm = diag(H);
-eta = abs(x-x_0in);
-h = x(2)-x(1);
-tot = m_order+s_order;
-S = [];
-M = [];
-poss = find(tot*h/2 >= eta);
+% Return zeros if x0 is outside grid
+if(x_0in < x(1) || x_0in > x(end) )
+
+    ret = zeros(size(x));
+
+else
+
+    fnorm = diag(H);
+    eta = abs(x-x_0in);
+    h = x(2)-x(1);
+    tot = m_order+s_order;
+    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);
+    elseif length(poss) == (tot + 1)
+        poss = poss(1:end-1);
+    end
+
+    % Use first tot grid points
+    if length(poss)<tot && eta(end)>eta(1)
+        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;
 
-% Ensure that poss is not too long
-if length(poss) == (tot + 2)
-    poss = poss(2:end-1);
-elseif length(poss) == (tot + 1)
-    poss = poss(1:end-1);
+    % Use last tot grid points
+    elseif length(poss)<tot && eta(end)<eta(1)
+        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
+    else    
+        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;
+    end
+
+    h_pol = pol(2)-pol(1);
+    b = zeros(m_order+s_order,1);
+
+    for i = 1:m_order
+        b(i,1) = x_0^(i-1);
+    end
+
+    for i = 1:(m_order+s_order)
+        for j = 1:m_order
+            M(j,i) = pol(i)^(j-1)*h_pol*norm(i);
+        end
+    end
+
+    for i = 1:(m_order+s_order)
+        for j = 1:s_order
+            S(j,i) = (-1)^(i-1)*pol(i)^(j-1);
+        end
+    end
+
+    A = [M;S];
+
+    d = A\b;
+    ret = x*0;
+    ret(index) = d/h*h_pol;
 end
 
-% Use first tot grid points
-if length(poss)<tot && eta(end)>eta(1)
-    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)
-    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
-else    
-    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;
-end
-
-h_pol = pol(2)-pol(1);
-b = zeros(m_order+s_order,1);
-
-for i = 1:m_order
-    b(i,1) = x_0^(i-1);
-end
-
-for i = 1:(m_order+s_order)
-    for j = 1:m_order
-        M(j,i) = pol(i)^(j-1)*h_pol*norm(i);
-    end
-end
-
-for i = 1:(m_order+s_order)
-    for j = 1:s_order
-        S(j,i) = (-1)^(i-1)*pol(i)^(j-1);
-    end
-end
-
-A = [M;S];
-
-d = A\b;
-ret = x*0;
-ret(index) = d/h*h_pol;
 end