t=linspace(0,2*pi,100);
f = [ones(1,50) -ones(1,50)];
f_3 = (4/pi)*sin(t) + (4/(3*pi))*sin(3*t);

f_21=zeros(1,100);
for k=1:2:21
    f_21 = f_21 + (4/(k*pi))*sin(k*t);
end

plot(t,f,t,f_3,t,f_21)
legend('f','f_3','f_{21}')