The function newton_interpolation takes in two arrays x and y as inputs, which represent the x and y coordinates of a set of data points. It then calculates the Newton Interpolation Polynomial of these data points, and displays both the polynomial equation and the divided difference table. The polynomial equation is displayed as a character stri…

newton_interpolation.m

function p = newton_interpolation(x, y)
    % Plots the Newton Interpolation Polynomial and the divided difference table.
    %   x - the x-coordinates of the data points
    %   y - the y-coordinates of the data points
    
    % Number of data points
    n = length(x);
    
    % Divided difference table
    F = divided_differences(x, y);
    
    % Newton polynomial    
    p = sym(F(1,1));
    for i = 2:n
        term = F(i,i) * prod(sym('x') - x(1:i-1));
        p = p + term;
    end
    disp(['Polynomial: ', char(p)]);
    
    % Return the polynomial as a function
    p = matlabFunction(p, 'Vars', {'x'});
    
    % Plot the results
    x_interp = linspace(min(x), max(x), 1000);
    y_interp = double(subs(p, 'x', x_interp));
    
    plot(x_interp, y_interp, 'b-', 'LineWidth', 1.5, 'Color', 'red');
    hold on;
    plot(x, y, 'o', 'MarkerFaceColor', 'red', 'MarkerEdgeColor', 'black');
    hold off;
    legend('Newton Interpolation Polynomial', 'Data Points');
    xlabel('X');
    ylabel('Y');
    title('NEWTON INTERPOLATION');
    grid on; box on; grid minor;
    set(gca,'FontSize',12,'LineWidth',1.5);
    
    % Display the divided difference table
    disp('Divided Difference Table:');
    disp(F);
end

function F = divided_differences(x, y)
    % Returns the divided difference table.
    %   x - the x-coordinates of the data points
    %   y - the y-coordinates of the data points

    % Number of data points
    n = length(x);
    F = zeros(n, n);
    F(:, 1) = y';
    for j = 2:n
        for i = j:n
            F(i, j) = (F(i, j-1) - F(i-1, j-1)) / (x(i) - x(i-j+1));
        end
    end
end