I agree that my comment on what the Gibbs phenomenon actually is was not very helpful to the actual problem at hand. The expression for the sawtooth is
A*( 1/2 - Sum((1/(pi*n))*sin(n*w0*t)) ).
where A is the double amplitude of the triangle wave, A = 4 in this case. Since the coefficients are A / (pi*n), your code is too low by a factor of 2. With the method you used, if you had included both positive and negative frequencies in the sum instead of just the positive ones, you would have got the factor of 2.
You have c0 = 2 in the code but it is not being added in successfully. Also your for loops are not exactly for loops since each set of instructions is only done once.
Making a new variable name for every value of n (dynamic variable naming) as you are doing is highly unrecommended. For one thing, if you want to do this process for, say, 100 terms, it’s just not practical. Rather than naming row vectors c1,c2 etc. it’s much better to put all the rows into a matrix and address each row with a row index. So you can do the following.
t = -10:.1:30;
c0 = 2;
y = zeros(N,length(t));
for n = 1:N
cn1 = 4j/(n*pi);
y(n,:) = real(cn1*exp(1j*n*wo*t));
yT = c0 + sum(y);
Now to do 50 terms you just change N to 50. I hope you will agree that this is highly efficient compared to making a new variable name for every term.
If you do change N to 50 you can see the actual Gibbs phenomenon which are the spikes at y = 4.3 and y = –0.3, approximately. Those get narrower with large N, but they don’t go away. Everywhere else the fourier sum gets closer and closer to the sawtooth wave.