As Thorny said, Milnor's axiomatic definition seems to be precisely the best way of proving that different definitions are the same. The main thrust of his "definition" is the proof that any invariants that satisfy these axioms must be the same as Stiefel-Whitney classes. In his book, they connect the two notions I describe below as well as the Steenrod-squares definition. They should also serve to prove that all the definitions you talk about are the same.
The rest of this answer might have less to do with your exact question than with my tendency to see an interesting question title and start writing. Sorry! Still, I feel that they are things that should be said (or, at least, don't deserve to be deleted).
I think there are two very important ways to understand characteristic classes. Both are explained in Milnor's Characteristic Classes, but not as the definition, since they are not as precise (but, to me, they are much more intuitive).
Think of your vector bundle as a map from your space X into a Grassmanian. The cohomology of the Grassmanian (more precisely, either the $\mathbb Z/2$ cohomology of the real Grassmanian, or the usual cohomology of the complex Grassmanian) is a polynomial algebra on some generators. The characteristic classes (Stiefel-Whitney or Chern, respectively) are precisely the pullbacks of these cohomology classes to X via the map.
Reading your question carefully, I guess you already knew this. Still, I think you should give this definition more credit. In particular, I think that this is the best explanation of the philosophical reason why "characteristic classes" exist. On thing that confuses me: why are the pullbacks of the integer cohomology of the real Grassmanian never called characteristic classes? I'm sure they are a pain to calculate, but that doesn't justify why nobody seems to care for them at all...
You can understand them through obstruction theory (another reference: Steenrod's "Theory of Fibre Bundles). The idea is to generalize the definition of the Euler characteristic using vector fields. Namely, try to construct a nowhere-zero section of your bundle. The obstruction will be a cohomology class, which is called the Euler class (and corresponds mod 2 to the top Stiefel-Whitney class). Try to construct two linearly independent nowhere-zero sections of the bundle. The obstruction will be a cohomology class which, mod 2, will be the next (one dimension lower) Stiefel-Whitney class. If you keep going like this, you'll construct all the classes.
Here is an explanation of why the obstructions to constructing non-zero sections are cohomology classes, for the case of a single section.
Think of your space X as a CW-complex; start constructing it on the 0-skeleton, and then try extending the section to 1-skeleton, and so on. At each step, you will basically be solving the following problem:
Given a vector field on the boundary $S^{n-1}$ of the ball $B^n$, can you extend it to the whole ball?
To solve this, think of the vector field as a map $S^{n-1}\to \mathbb R^m$ where m is the dimension of your bundle (you can assume that the bundle is trivial over the ball $B^n$ since the ball is contractible). Since the vector field is supposed to be nowhere zero, you can think of this as a map $S^{n-1}\to S^{m-1}$. If $n<m$, this map is always nullhomotopic and always extends to the ball. If $n=m$, you get an integer, the degree of the map, which tells you if you can extend. Since you get an integer for each degree-m cell of the CW-complex, you get something that looks like a cohomology class in $H^m(X)$ (of course, you need to verify separately that it actually is one, and if you are precise enough, you'll see that these integers only make sense mod 2). This is the Euler class.
If you wanted to construct two linearly independent sections, first construct one up to the $n-1$-skeleton (which is always possible). Now, let's start making the second one. You might as well require the second section to be orthogonal to the first. So, in the extension problem, you'll have a map $S^{n-1}\to \mathbb R^{m-1}$ where the $\mathbb R^{m-1} \subset \mathbb R^m$ is the subspace orthogonal to the first section. Since it also can't be zero, it's really a map $S^{n-1}\to S^{m-2}$. The rest of the argument is the same; you get a class in $H^{m-1}(X)$.
Usual disclaimer: there may be mistakes anywhere. Please point them out!
I'm grateful to Allen Hatcher, who pointed out that this answer was incorrect. My apologies to readers and upvoters. I thought it more helpful to correct it than delete outright, but read critically.
If $X$ and $Y$ are cell complexes, finite in each degree, and two maps $f_0$ and $f_1\colon X\to Y$ induce the same map on cohomology with coefficients in $\mathbb{Q}$ and in $\mathbb{Z}/(p^l)$ for all primes $p$ and natural numbers $l$, then they induce the same map on cohomology with $\mathbb{Z}$ coefficients. To see this, write $H^n(Y;\mathbb{Z})$ as a direct sum of $\mathbb{Z}^{r}$ and various primary summands $\mathbb{Z}/(p^k)$, and note that the summand $\mathbb{Z}/(p^k)$ restricts injectively to the mod $p^l$ cohomology when $l\geq k$. One can take only those $p^l$ such that there is $p^l$-torsion in $H^\ast(Y;\mathbb{Z})$. (I previously claimed that one could take $l=1$, which on reflection is pretty implausible, and is indeed wrong.)
We can try to apply this to $Y=BG$, for $G$ a compact Lie group. For example, $H^{\ast}(BU(n))$ is torsion-free (and Chern classes generate the integer cohomology), and so rational characteristic classes suffice. In $H^{\ast}(BO(n))$ and $H^{\ast}(BSO(n))$ there's only 2-primary torsion. That leaves the possibility that the mod 4 cohomology contains sharper information than the mod 2 cohomology. It does not, because, as Allen Hatcher has pointed out
in this recent answer,
all the torsion is actually 2-torsion.
It's sometimes worthwhile to consider the integral Stiefel-Whitney classes $W_{i+1}=\beta_2(w_i)\in H^{i+1}(X;\mathbb{Z})$, the Bockstein images of the usual ones. These classes are 2-torsion, and measure the obstruction to lifting $w_i$ to an integer class. For instance, an oriented vector bundle has a $\mathrm{Spin}^c$-structure iff $W_3=0$.
[I'm sceptical of your example in $2\mathbb{CP}^2$. So far as I can see, $3a+3b$ squares to 18, not 6, and indeed, $p_1$ is not a square.]
Best Answer
Here is one way I like to think of the analogy.
The maximal torus of diagonal matrices $T^{n} \subset U(n)$ gives a map $BT^n \to BU(n)$ which on integral cohomology gives an isomorphism from $H^{\ast}(BU(n))$, which is a polynomial algebra on $n$ generators of degrees $2i$, $1\leq i\leq n$, to $H{\ast}(BT^n)^{\Sigma_n}$, which is the polynomial algebra on the symmetric polynomials $\sigma_i$ in the $n$ standard degree 2 generators. The Chern class $c_i$ is the element of the domain that maps to $\sigma_i$.
The maximal $2$-torus of diagonal matrices $(C_2)^{n} \subset O(n)$ gives a map $BC_2^{n} \to BO(n)$ which on mod $2$ cohomology gives an isomorphism from $H^{\ast}(BO(n))$, which is a polynomial algebra on $n$ generators of degrees $i$, $1\leq i\leq n$, to $H^{\ast}(BC_2^n)^{\Sigma_n}$, which is the polynomial algebra on the symmetric polynomials $\sigma_i$ in the $n$ standard degree $1$ generators. The Stiefel-Whitney class $w_i$ is the element of the domain that maps to $\sigma_i$.
Thinking of $BT^{n}$ as $(CP^{\infty})^{n}$ and $BC_2^{n}$ as $(RP^{\infty})^{n}$ may help.