Plots the Newton polynomial interpolating a given function (and shows the Runge phenomenon)

plot_newton.m

function plot_newton
f = @(x) 1./(1+x.^2);   % function to interpolate
x = -5:5;               % nodes of interpolation
y = f(x);
c = coeff(x,y);         % calculate Newton polynomial coefficients
t = linspace(-5,5);
poly = newton(c,x,t);   % evaluate the polynomial 
plot(t,[f(t);poly]);    % plot function and its polynomial
hold on 
plot(x,y,'ro')          % also mark the interpolated values
end

function c = coeff(x,y)
n = length(x);
for k = 2:n
    y(k:n) = (y(k:n)-y(k-1))./(x(k:n)-x(k-1));
end
c = y;
end

function p = newton(c,x,t)
n = length(x);
p = c(n)*ones(size(t));
for k = n-1:-1:1
    p = (t-x(k)).*p + c(k);
end
end