% Continuous least squares approximation of sqrt(x) by polynomials for x in [0,1].
% Construction uses the monomial basis, giving a numerically treacherous Hilbert matrix.

 N = input('polynomial degree? ');

 A = zeros(N+1,N+1);
 f = zeros(N+1,1);

 for j=0:N, 
    for k=0:N
       A(j+1,k+1) = 1/(j+k+1);
    end
    f(j+1) = 1/(j+1.5);         % integral of sqrt(x)*x^j from 0 to 1 
 end

 p = A\f; 

 figure(1), clf
 xx = linspace(0,1,500);
 plot(xx,sqrt(xx), 'b-','linewidth',2), hold on
 plot(xx,polyval(flipud(p),xx), 'k--','linewidth',2)
 legend('f(x) = sqrt(x)', 'poly approx',2)
 xlabel('x','fontsize',20)
 fprintf('Polynomial coefficients, monomial basis\n')
 disp([p])