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

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

% Use last tot grid points
elseif length(poss)<tot && eta(end)<eta(1)
    index = length(x)-tot:tot;
    pol = (x(end-tot:end)-x(end-tot))/(x(end)-x(end-tot));
    norm = fnorm(end-tot:end)/h;
    x_0 = (x_0in-x(end-tot))/(x(end)-x(end-tot));

% 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