Polynomial interpolation using Chebyshev nodes

chebyshev_interpolation.m

function chebyshev_interpolation
a = -5;
b = 5;
t = linspace(a,b);
plot(t, [chebyshev(a,b,5,t); chebyshev(a,b,10,t); chebyshev(a,b,15,t); chebyshev(a,b,30,t)]);
hold on
plot(t, f(t), 'k--');  % original function in dashed black curve
legend('degree 5', 'degree 10', 'degree 15', 'degree 30', 'function')  % describe the curves
end

function y = f(x)  % using named function, so that it can be accessed from elsewhere
y = 1./(1+x.^2);   % an anonymous function would be local to some function block  
end

function p = chebyshev(a,b,n,t)
k = 1:n;
x = (a+b)/2+(b-a)/2*cos(pi*(2*k-1)/(2*n)); % all nodes computed at once
c = coeff(x, f(x));
p = newton(c, x, t);
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