% written by Mary Pugh (mpugh@math.toronto.edu) on Oct 3, 2002.
clear
% choose the interval
x = 0:.01:2*pi;

% choose a function
% f = x.*(2*pi-x);
f = sin(8*x);
% f = exp(x);

% random f:
a = rand(1,20);
b = 2*pi*rand(1,20);
% f = zeros(size(x));
% for k=1:20;
%   f = f + a(k)*cos(k*x + b(k));
% end

% make the range of f be in [0,1]:
f = f - min(f);
f = f/max(f);

% save the original function
F = f;
% the number of approximants
N = 10;

for n=1:N
  g(n,:) = max(0,min(f-2^(n-1)/3^n,2^(n-1)/3^n));
  f = f - g(n,:);
end

figure(1);
clf
title('f is in yellow, the partial sums are in red')
plot(x,g(1,:),'r')
for n=2:N
  plot(x,F);
  hold on;
  plot(x,sum(g(1:n,:)),'r')
  hold off;
  pause(2)
end


figure(2);
title('here are the g_n functions')
for n=1:N
  plot(x,g(n,:))
  pause(2)
end


disp('here are the errors')
format compact
format short e
err(1) = max(abs(F-g(1,:)));
for n=2:N
  err(n) = max(abs(F-sum(g(1:n,:))));
end
err
disp('they are decreasing by 2/3 at each step')