understanding Steenrod squares
Here's one way to understand them. The external cup square $a \otimes a \in H^{2n}(X \times X)$ of $a \in H^n(X)$ induces a map $f:X \times X \to K(Z_2, 2n)$. It can be show that this map factors through a map $g:(X \times X) \times_{Z_2} EZ_2 \to K(2n)$, where $Z_2$ acts on the product by permuting the factors and $EZ_2$ can be taken to just be $S^\infty$. If you unravel what this means, it says that our original map $f$ was homotopic to the map obtained by first switching the coordinates and then applying $f$. It also says that this homotopy, when applied twice to get a homotopy from $f$ to itself, is homotopic to the identity homotopy, and we similarly have a whole series higher "coherence" homotopies. Now $X \times BZ_2$ maps to $(X \times X) \times_{Z_2} EZ_2$ as the diagonal, so we get a map $X \times BZ2 \to K(2n)$. But $BZ_2$'s cohomology is just $Z_2[t]$, so this gives a cohomology class $Sq(a) \in H^*(X)[t]$ of degree $2n$. If we write $Sq(a)=\sum s(i) t^i$, it can be shown that $s(i)=Sq^{n-i}a$.
What does this mean? Well, if our map $f$ actually was invariant under switching the factors (which you might think it ought to be, given that it appears to be defined symmetrically in the two factors), we could take $g$ to be just the projection onto $X \times X$ followed by $f$. This would mean that $Sq(a)$ comes from just projecting away the $BZ_2$ and then using $a^2$, i.e. $Sq^n(a)=a^2$ and $Sq^i(a)=0$ for all other $i$. Thus the nonvanishing of the lower Steenrod squares somehow measures how the cup product, while homotopy-commutative (in terms of the induced maps to Eilenberg-MacLane spaces), cannot be straightened to be actually commutative. Indeed, in the universal example $X=K(Z_2,n)$, the map $f$ is exactly the universal map representing the cup product of two cohomology classes of degree $n$.
Some somewhat terse notes on this can be found here; see particularly part III. (Sorry, the link is now dead.)
Here's how I explain Steenrod squares to geometers. First, if $X$ is a manifold of dimension $d$ then one can produce classes in $H^n(X)$ by proper maps $f: V \to X$ where $V$ is a manifold of dimension $d-n$ through many possible formalisms - eg. intersection theory (the value on a transverse $i$-cycle is the count of intersection points), or using the fundamental class in locally finite homology and duality, or Thom classes, or as the pushforward $ f_*(1) $ where $1$ is the unit class in $H^0(V)$. Taking this last approach, suppose $f$ is an immersion and thus has a normal bundle $\nu$. If $x = f_*(1) \in H^n(X)$ then $Sq^i(x) = f_*(w_i(\nu))$. This is essentially the Wu formula.
That is, if cohomology classes are represented by submanifolds, and for example cup product reflects intersection data, then Steenrod squares remember normal bundle data.
You can understand the squares purely algebraically: let $R = F_2[x_1, x_2, \ldots, x_n]$ be a polynomial ring over the field $F_2$ of $2$ elements. A ring homomorphism $f$ on $R$ is completely determined by the values of $f(x_i)$, and those values are unrestricted, right? So define $f \colon R \to R$ by $f(x) = x + x^2$ for each $x = x_i$. (Since $\operatorname{char}(F)=2$, we actually have $f(x) = x + x^2$ for each $F$-linear combination of the $x_i$, too, i.e. this is a basis-free -- "gentlemanly"? -- definition.) Then for any homogeneous polynomial $p$ in $R$, we have $f(p) = p + p^2 + \text{other stuff}$; separating the components of $f(p)$ by degree, we can write $f(p) = \sum Sq^i ( p )$. E.g. $f(x_1 x_2) = (x_1 + x_1^2) \cdot (x_2 + x_2^2) = (x_1 x_2) + (x_1^2 x_1 + x_2^2 x_1) + (x_1 x_2)^2$, so $Sq^1( x_1 x_2 ) = x_1 x_2 (x_1 + x_2)$.
You can extend the definition of the $Sq^i$ to non-homogeneous polynomials by additivity if you like. This definition of the $Sq^i$ is consistent for polynomials in arbitrary numbers of variables (i.e. there are commutative diagrams involving the inclusions $F_2[x_1] \to F_2[x_1,x_2]$ etc.) So you can visualize the whole Steenrod algebra as an algebra of endomorphisms on an infinite polynomial ring.
This definition is sufficient to check the Adem relations.
Topologically, this mechanism defines $Sq^i$ on $H^\ast( (\mathbb{RP}^\infty)^n, Z/2Z)$. Then by naturality it defines the $Sq^i$ on cohomology classes that are pullbacks from those rings, i.e. you have a description of how $Sq^i$ acts on cohomology classes defined by vector bundles. This really is how I think of them when I play with classifying spaces $BG$.