 x = chebfun('x',[0 1]);
 f = exp(x).*sin(5*pi*x);   % the function to approximate
 f = exp(x);

 n = 4;                     % the degree of the interpolating polynomial
 X = [0:n]/n;               % interpolation points

 clear L
 figure(1), clf
 plot([0 1], [0 0], 'r-','linewidth',2)
 hold on
 for k=1:n+1
    plot(X(k)*[1 1], [-.5 1.5],'-','color',.75*[1 1 1])
 end

 for k=1:n+1
    L{k} = chebfun('1',[0 1]);
    for i=1:n+1
       if i~=k, L{k} = L{k}.*((x-X(i))/(X(k)-X(i))); end
    end
    plot(L{k},'k-','linewidth',2)
    hold on 
%    pause
 end


 pn = chebfun('0',[0 1]);
 for j=1:n+1
    pn = pn + f(X(j))*L{j};
 end

 figure(2), clf
 plot(f,'b-','linewidth',3) 
 hold on
 plot(pn,'r--','linewidth',3) 
 plot(X,f(X),'bo','linewidth',2,'markersize',14) 
