% Math 505; Gautchi 2.24: Lebesgue functions
fprintf('        n         equispaced            Chebyshev \n')
for ii=0:2
    n=5*2^ii; xa(1:n+1,1)=0; xb(1:n+1,1)=0;
   for i= 1:n+1
      xa(i,1) = -1+2*(i-1)/n;
      xb(n+2-i,1) = cos(pi*(2*(i-1)+1)/(2*n+2));
   end
   for i=1:n
      for j=1:20
         ga(20*(i-1)+j,1)=((21-j)*xa(i)+(j-1)*xa(i+1))/20;
         gb(20*(i-1)+j,1)=((21-j)*xb(i)+(j-1)*xb(i+1))/20;
      end
   end
   ga(20*n+1,1)=xa(n+1);
   gb(20*n+1,1)=xb(n+1);
   for k = 1:20*n+1                    %compute Lebesgue functions
      na(:,k)=ga(k)-xa(:);    nb(:,k)=gb(k)-xb(:);       %numerators
   end
   for i = 1:n+1                           % denominators
       for j=1:n+1
           da(i,j) = xa(i)-xa(j);   db(i,j) = xb(i)-xb(j);
       end
   end
for k = 1:20*n+1
      for i = 1:n+1
          la(i,k) = 1; lb(i,k) = 1;
          for j=1:n+1
              if j~=i
                  la(i,k)=la(i,k)*na(j,k)/da(i,j);   lb(i,k)=lb(i,k)*nb(j,k)/db(i,j);
              end
          end
      end
      lama(k)=1;  lamb(k)=1;
      for i=1:n+1
          lama(k)=lama(k)+abs(la(i,k));       lamb(k)=lamb(k)+abs(lb(i,k));
      end
   end
   figure
   subplot(2,1,1);    plot(ga,log(lama))
   subplot(2,1,2);    plot(gb,lamb)
   lamax =lama(1);   lambx =lamb(1); % compute Lebesgue constants
   for k = 2:20*n+1
       if lama(k) > lamax;   lamax= lama(k);    end
       if lamb(k) > lambx;   lambx= lamb(k);    end
   end 
   fprintf('%9i %22.14f %22.14f \n',n,lamax,lambx);  clear
end
