function interpolate

% lagrange interpolation of f(x) for a < x < b
% f(x) defined at end of file
% uses nx data points

nx=5
a=-pi; b=pi;

% compute data points
xd=linspace(a,b,nx);
for ix=1:nx
    yd(ix)=f(xd(ix));
end

% evaluate interpolation poly at n points in interval
n=1000;
x=linspace(a,b,n);
for ii=1:n
    p(ii)=0;
    for k=1:nx
        p(ii)=p(ii)+yd(k)*ell(k,x(ii),xd);
    end
end

% compute function
for ie=1:n
    ye(ie)=f(x(ie));
end

% plot results
clf
% get(gcf)
set(gcf,'Position', [1 925 560 420])
hold on
box on
plot(x,ye,'r','LineWidth',1.5)
plot(x,p,'--b','LineWidth',1.5)
xlim([a b])
say=['Number of Data Points = ',num2str(nx)];
title(say,'FontSize',20,'FontWeight','bold')
grid on
xlabel('x-axis')
ylabel('y-axis')
legend({' Exact',' Lagrange'},'Location','South','FontSize',14,'FontWeight','bold')
set(gca,'FontSize',14,'FontWeight','bold')
hold off

% lagrange basis function
function p=ell(i,x,xd)
[n1 n2]=size(xd);
p=1;
for j=1:n2
    if j ~= i
        p=p*(x-xd(j))/(xd(i)-xd(j));
    end
end

% the function f(x)
function g=f(x)
g=sin(x);