%-----------------------------------------------------------------------
% Find polynomial approximation for a given function or data points
%-----------------------------------------------------------------------
% Instructor: Nam Sun Wang
%-----------------------------------------------------------------------

% Start fresh ----------------------------------------------------------
      clear all

% Program header -------------------------------------------------------
      disp('This program finds the interpolation polynomial to the following function:')
      disp('  f=exp(-x^2)')

% Generate interpolation points ----------------------------------------
% n points (i.e., n-1 intervals) between x=[alpha,beta]
      n=input('Enter number of interpolating points: ');
      if n<=0 n=11; end
      Chebychev = input('Use Chebychev points (y/N): ', 's');

      alpha = -5;
      beta  =  5;
      x=zeros(n,1);
      y=zeros(n,1);
      for i=1:n
        if Chebychev=='y' | Chebychev=='Y'
          x(i) = (alpha+beta)/2 + (alpha-beta)/2*cos((i-1)/(n-1)*pi);   %Chebychev points
        else
          x(i) = alpha + (beta-alpha)/(n-1)*(i-1);   %Equi-distant points
        end
        y(i) = exp(-x(i)^2);
      end

% Find the interpolation polynomial coefficients -----------------------
      for j=1: n
        Z(:,j) = x.^(j-1);
      end
      a = Z^(-1)*y;

% Plot the given function and the interpolation polynomial -------------
      xx=[alpha: 0.05: beta]';
      yy=zeros(size(xx));
      for i=1:size(xx,1)
        f(i) = exp(-xx(i)^2);
        yy(i)=a(1);
        for j=2:n
          yy(i) = yy(i) + a(j)*xx(i)^(j-1);
        end
      end
      plot(xx, f, '-', xx, yy, '-', x, y, 'o')
      title('Interpolation Polynomial for f=exp(-x^2)')
      legend('Given function', 'Interpolating Polynomial', 'Interpolation Points')