function [coefs,x,P] = ChebyshevPol(f,a,b,n,varargin)
% ChebyshevPol('f',a,b,n) returns the vector of coefficients
% of the Chebyshev Polynomial approximation of f in the interval [a,b].
% n is the desired order of the polynomial approximation. 
%
% [c, x, P] = ChebyshevPol(f,a,b,n) returns the coefficients, the
% Gauss-Chebyshev nodes x and a matrix P with the i-th
% Chebyshev polynomial at values x in column i+1. 
%
% WARNING: coefficients refer to Chebyshev Polynomials over [-1, 1]. 
% To evaluate over x in [a,b], transform to z = 2*(x-a)/(b-a)-1 before
% you apply ChebyshevPol.

% copyright 2001 Leonardo Rezende.
% This program can be freely copied, used and modified, as long as this
% copyright notice is kept in place. Please inform the author of any
% useful modifications, at lrezende@stanford.edu.

m = n+1;

x = cos(pi*((1:m)-1/2)/m)';
w = pi/n*ones(m,1);

% calculates Chebyshev Polynomials

P = [ones(size(x)), x];
for j = 2:n;
    p = 2*x.*P(:,j)-P(:,j-1);
    P = [P, p];
end;

% evaluates the function at nodes

F = [];
z = (b-a)/2*(x+1)+a;
for j = 1:size(x)
  F = [F; feval(f, z(j), varargin{:})];
end

coefs = P'*diag(w)*F ./ diag(P'*diag(w)*P);

