Try this:
f = 5; % Frequency
N = 16; % Order Of Harmonics
t = linspace(0, 4, 500); % Time Vector
sqwv = @(f,t,n) 4*sum(bsxfun(@rdivide, sin(2*pi*(1:2:N)'*t), (1:2:N)'))/pi;
figure
plot(t, sqwv(f,t,N))
grid
