Rotational oblateness I am trying to compute the amount of oblateness that is caused by planetary rotation. I picture the force of gravity added to the centrifugal force caused by the rotation of the planet as follows:
$\hspace{4cm}$
That is, at the point in question, at latitude $\phi$, the distance from the axis of rotation is $r\cos(\phi)$. Thus, the centrifugal force would be $\omega^2r\cos(\phi)$ in a direction perpendicular to the axis of rotation. The radial and tangential components would be $\omega^2r\cos^2(\phi)$ and $\omega^2r\cos(\phi)\sin(\phi)$, respectively.
My assumption is that the surface of the planet would adjust so that it would be perpendicular to the effective $g$; that is, the sum of the gravitational and centrifugal forces. This would lead to the equation
$$
\frac{\mathrm{d}r}{r\,\mathrm{d}\phi}=-\frac{\omega^2r\cos(\phi)\sin(\phi)}{g-\omega^2r\cos^2(\phi)}.
$$
We can make several assumptions here, and I will assume that $\omega^2r$ is small compared to $g$. Thus, we get
$$
\int_{\text{eq}}^{\text{np}}\frac{\mathrm{d}r}{r^2}
=-\frac{\omega^2}{g}\int_0^{\pi/2}\cos(\phi)\sin(\phi)\,\mathrm{d}\phi
$$
which leads to
$$
\frac1{r_{\text{np}}}-\frac1{r_{\text{eq}}}
=\frac{\omega^2}{2g}
$$
and
$$
1-\frac{r_{\text{np}}}{r_{\text{eq}}}
=\frac{\omega^2r_{\text{np}}}{2g}.
$$
However, numerical evaluation and Wikipedia seem to indicate that this should be twice what I am getting. That is,
$$
1-\frac{r_{\text{np}}}{r_{\text{eq}}}
=\frac{\omega^2r^3}{Gm}
=\frac{\omega^2r}{g}.
$$
What am I doing wrong?
 A: As is well-known from Newton's shell theorem, the gravitational field $g(r)=\frac{GM}{r^2}$ outside a spherically symmetric mass-distribution is the same as if the total mass $M$ sat in the center.
It seems that OP wants to calculate the oblateness of Earth under the simplifying assumption that the backreaction (which the re-distributed mass has on Earth's gravitational field) can be ignored. In other words, we assume that the gravitational field is given by just the monopole contribution $g(r)=\frac{GM}{r^2}$, and we neglect higher multipole moments in a multipole expansion.
I) This is what Mark Eichenlaub did in this Phys.SE post. To compare let us replace the latitude $\phi$ with the polar angle $\theta=\frac{\pi}{2}-\phi$. The total potential energy is a sum of the gravitational monopole potential energy and the centrifugal potential energy
$$ U~=~- \frac{GM}{r}-\frac{(\omega r \sin \theta)^2}{2}. \tag{1}$$
The point is now that (in this idealized model) the surface of Earth is an equipotential surface. ("Else the water in the oceans would rush to re-distribute itself.")
Comparing Equator and the North pole leads to
$$ g(a)h~\approx~\frac{GM}{b}-\frac{GM}{a}~\stackrel{(1)}{=}~\frac{(\omega a)^2}{2}~>~0,\tag{2} $$
where $a$ and $b$ are the equatorial and polar radius of the Earth, respectively; and $h:=a-b\ll a$ is the sought-for height difference. Equation (2) leads precisely to Mark Eichenlaub's monopole result for $h$, which is $\frac{2}{5}$ smaller than the quadrupole result.
II) On the other hand, if we differentiate eq. (1), we get precisely OP's force equilibrium formula
$$ 0~=~\mathrm{d}U~=~\left(g(r)-(\omega\sin \theta)^2 r\right)\mathrm{d}r -(\omega r )^2\sin \theta \cos \theta \mathrm{d}\theta. \tag{3}$$
At this point OP ignores the radial dependence of $g(r)$, and treats it as a constant $g$. This model corresponds to a total potential energy
$$ V~=~ gr-\frac{(\omega r \sin \theta)^2}{2}. \tag{4}$$
Comparing Equator and the North pole leads to
$$ gh~=~g b-ga~\stackrel{(4)}{=}~-\frac{(\omega a)^2}{2}~<~0, \tag{5}$$
which predicts an prolate Earth rather than an oblate Earth.
III) Next OP assumes that one of the centrifugal terms $(\omega\sin \theta)^2 r \ll g$ in eq. (3) is small and should be ignored. This mean that eq. (3) is no longer a perfect differential. However an integrating factor is $\frac{1}{r^2}$, so this leads to a first integral
$$  W~=~ -\frac{g}{r}-\frac{(\omega \sin \theta)^2}{2}. \tag{6}$$
Comparing Equator and the North pole leads to
$$ \frac{gh}{ab}~=~\frac{g}{b}-\frac{g}{a}~\stackrel{(6)}{=}~\frac{\omega^2}{2}~>~0, \tag{7} $$
which remarkably reproduces Mark Eichenlaub's monopole result (2). In other words, two not-so-small-approximations by OP have cancelled out.
A: Problem is with the assumption:

... I will assume that $\omega^{2}r$ is small compared to $g$ 

This means: centrifugal force at the equator is negligible.
And this cannot be true, because the planet would be perfect sphere.
(Note: word "small" was used as "negligible" when constructing the equation)
