c_0 = 1;
N = 40;

% define the coefficients of the powerseries for sqrt(1-t) about t=0:
c(1) = -1/2;
c(2) = -1/8;
for n=3:N
  c(n) = -(gamma(2*n-2+1)/(2^(n-1)*gamma((n-1)+1)*gamma(n+1)*2^n));
end

s = -1:.01:1;
s = s';
% choose the function whose absolute value you want to approximate
f = s;
% choose the function whose absolute value you want to approximate
f = sin(5*s);
%n = length(s);
%g = rand(n+2,1);
% choose a random function and smooth it a little...
%f = (g(1:n) + g(2:n+1) + g(3:n+2))/3;
%f = f*(f-1/2);

% make sure the L^infty norm of f is 1.
M = max(abs(f));
f = f/M;


t = 1-f.^2;
% define the approximants as in the proof of Stone-Weierstrauss
p(:,1) = c_0*ones(size(t));
for n=2:N+1
  p(:,n) = p(:,n-1) + c(n-1)*t.^(n-1);
end

disp('number of approximants: '),N

clf;
figure(1);
for n=1:N+1
  plot(s,abs(f),'g');
  hold on
  plot(s,p(:,n),'r');
  hold off;
  title('in green: |f|, in red: p_n')
  figure(1);
  pause(1);
end

clf;
figure(2);
for n=1:N+1
  plot(s,p(:,n)-abs(f))
  title('plot of |p_n(x) - f(x)|')
  figure(2);
  pause(1);
end