What is the formula for a square wave? After writing this, I realised it's rather long. Please delete if too convoluted.
Many many years ago I studied electronics and in a class we used excel to plot a sine wave.
Simple.
Get a sample rate (24) and make a column of 0 to 360 in 15° increments.
Next column, 360 divided by the previous column (decimal point form in 1/24 increments) and multiply by 2 (for PI)
Next column, multiply previous column by PI.
Next column, apply SIN() to previous column.
Here, we can apply a line chart, smooth it out and it looks like a reasonable sine wave.

Then my teacher told me at the time that a square wave is just a sine wave with harmonics added to it, specifically every odd harmonic.
Now onto making the harmonic column. A 3rd of the amplitude and 3 times the frequency, 5th and so on...
And we get a graph like this (albeit, rough):

Then we add all of them up at each 15° interval:

Now, it's sort of looking like a square wave but I was expecting it to be more... Square?
I then decided to revisit it and did it all the way up to the 21st harmonic and in 1° increments:
21st harmonic

Fundamental + 3rd - 21st (every odd) added

Is this true? Wikipedia says that it is, but I can't really see it?

The ideal square wave contains only components of odd-integer harmonic frequencies

 A: This is called the Wilbraham-Gibbs phenomenon: this "ringing" occurs for any discontinuous function (the artefacts you see in low-quality JPEGs are also related to this).
For functions that have a finite number of finite discontinuities, you can improve convergence considerably by using Fejér sums of the series: this causes the initial terms to dominate, so the oscillations from the terms with shorter wavelengths are relatively suppressed.
A: It seems you have two questions:  first, whether Fourier analyzing a square wave gives only odd harmonics and whether the approach you are following will converge to a square wave.  The answer to both is yes.  For the first, all the odd harmonics are odd functions of time, while the even harmonics are even functions of time.  If you Fourier analyze any odd function you will use only the odd harmonics.  The square wave is odd in time, so it uses only odd harmonics. For the second, your last plot is very useful.  It should look to you like it is converging nicely on a square wave in the middle.  Away from the transitions it is quite constant.  Those little ripples will disappear as more terms are added.  The worrisome thing is the big bumps near the transition.  As you add more terms they actually get taller and narrower and move towards the transition points.  For any given time, the peaks will eventually move between it and the transition and the value of the expansion will settle down to $\pm 1$ as you want.  
What is going on here is that sine and cosine functions are infinitely differentiable, so they don't cope well with discontinuities.  The coefficients in the Fourier expansion of a square wave fall off as $\frac 1n$, as they do for any discontinuous function.  If you expand a continuous function they will eventually decrease as $\frac 1{n^2}$.  If you expand a once differentiable function they will fall as $\frac 1{n^3}$  You can look at a table of Fourier expansions to see this.
