function [a0,An,Bn]=mysquare(T,N)
syms t n
w0=2*pi/T;
a0=1/T*(int(1,t,0,T/2));
an=2/T*(int(1*cos(n*w0*t),t,0,T/2));
bn=2/T*(int(1*sin(w0*n*t),t,0,T/2));
An=symsum(an.*cos(n*w0*t),n,1,N);
Bn=symsum(bn.*sin(n*w0*t),n,1,N);
%main script
T=4;
[a0,an,bn]=mysquare(T,10000);
ft=a0+an+bn;
fplot(ft);
This works fine for N<=1000 but for larger values it takes alot of time and no result.
Best Answer