function NDD = newton_divided_difference(x, y, a)

    n = length(x);

    diff_table = zeros(n, n);

 

    for i = 1 : n

        diff_table(i, 1) = y(i);

    end

 

    for j = 2 : n

        for i = 1 : n - j + 1

            diff_table(i, j) = (diff_table(i + 1, j - 1) - diff_table(i, j - 1)) / (x(i + j - 1) - x(i));

        end

    end

 

    NDD = diff_table(1, 1);

    val = (a - x(1));

 

    for j = 2 : n

        NDD = NDD + (val * diff_table(1, j));

        val = val * (a - x(2));

    end