%% script ChebyshevApproximation
%% compute the Chebyshev approximation of order n of the function fun
n = 2; 
accuracy = 1e-8;
A = 0; % left endpoint
B =+1; % right endpoint
fun = @(x)atan(x);

%%%%%%%%%%%%%% no modifications beyond this line %%%%%%%%%%%%
newFun = @(x)feval(fun,A+0.5*(x+1)*(B-A));
c = zeros(n+1,1);
k = 0;
c(1) = quad(@(p)feval(newFun,(cos(p)))*cos(k*p),0,pi,accuracy)*1/pi;
for k = 1:n
  c(k+1) = quad(@(p)feval(newFun,(cos(p)))*cos(k*p),0,pi,accuracy)*2/pi;
end%for

coeff = zeros(n+1,n+1);
for i=1:n+1;
  coeff(i,n-i+2:n+1) = Chebyshev(i-1);
end%for
newPol = coeff'*c

%%% solve for the coefficients on the interval [A,B]
y = polyval(newPol,linspace(-1,1,n+1))';
F = vander(linspace(A,B,n+1)');
% x = linspace(A,B,n+1)';
% F = zeros(n+1,n+1);
% for k = 0:n
%   F(:,k+1) = x.^(n-k);
% end%for
p = F\y


%%%%%%%%%%%%%%% visualize the results %%%%%%%%%%%%%%%%%%%%%
x  = linspace(A,B,101);
y1 = fun(x);
y2 = polyval(p,x);
figure(2);
plot(x,y2-y1)
grid on
title('absolute error of Chebyshev approximation')

figure(1);
plot(x,y1,x,y2)
legend('off')
grid on
title('function and Chebyshev approximation')


%%%%%% linear Regression, for comparison
F = zeros(length(x),n+1);
for k= 0:n
  F(:,k+1) = x'.^(n-k);
end%for

pLin = LinearRegression(F,y1');
y3 = polyval(pLin,x);
figure(3);
plot(x,y3-y1,x,y2-y1)
legend('linear regression','Chebyshev approximation')
%legend('show )
grid on
title('comparison of errors')

accuracyBit = -log2(max(abs(y2-y1))/max(abs(y1)))