a = -1;
b = 1;
x = linspace(a, b, 1000);
% f = @(x) sin(pi*x);
f = @(x) 1 ./ (1 + (5*x).^2);
n = 2;
xInterp = linspace(a, b, n+1);
p = interpolate(f, xInterp, x);

plot(x, p);
hold on;
plot(x, f(x), 'r');
plot(xInterp, f(xInterp), 'r*');
hold off;