% This program plots partial sums of the Fourier series of f(x)=1_{0,1/2}
% And averages of partial sums

x = (0:0.0001:1);
Lmax = 100;

Phi = zeros(Lmax+1,length(x));
SN = zeros(1,length(x));
sigmaN = zeros(1,length(x));
f = zeros(1,length(x));

f(x<0.5) = 1;

Phi(1,:) = 0.5 + (2/pi)*sin(2*pi*x);

for n=1:Lmax
    Phi(n+1,:) = Phi(n,:) + (2/(pi*(2*n+1)))*sin(2*pi*(2*n+1)*x);
end

SN = Phi(Lmax+1,:);
sigmaN = sum(Phi,1)/(Lmax+1);


figure(1)
plot(x,f,'g',x,SN,'r',x,sigmaN,'b')
grid on
legend('f','Symmetric partial sum: S_N','Cesaro mean: \sigma_N')
title(sprintf('Approximations for N=%d',Lmax+1))
