The first step is to find the partial derivative wrt r of the exterior gravitational potential of the oblate, derived by Gauss and Dirichlet: $$P = \frac{3 GMm}{(2 a e)} \left[B \left(1 - \frac{r^2}{(2 a^2 e^2)} + \frac{z^2}{(a^2 e^2}\right) + \frac{r^2}{(2 a^2 e^2)} \sin(B) \ \cos(B) - \frac{z^2}{(2 a^2 e^2)} \tan(B) \right] \ \ $$ where r and z are cylindrical coordinates, but I am ignoring the longitudinal $\phi$ azimuth coordinate due the symmetry. r and z are effectively cartesian coordinates of the elliptical cross section, where only positive values of r and z are valid.
The next step is to substitute the value of z (The vertical coordinate) as a function of r (the horizontal coordinate).
$$z(r) = \sqrt{(a^2 - a^2 e^2)*(1 - r^2/a^2)}$$
which is readily obtained from the definition of an ellipse.
A similar procedure to find the partial derivative wrt to z is used to find the vertical component Fz of the force:
The exact result for the horizontal component $F_r$ of the surface gravitational force can be seen in this Mathematica worksheet , but I recommend the simplified equations found by @Ghoster in his answer.
In short $F_z \propto M \sqrt{a^2-r^2}$ and $F_r \propto M r$.
While it may be initially surprising to see that the horizontal component (blue curve below) increases with distance from the centre, it becomes obvious this must be the case, since centrifugal force (CF) also increases linearly with increasing distance from the centre as per $CF \propto \omega^2 r$ and a fluid oblate rotating body is the reult of the equilibrium of those forces.
The resultant total force Ft is easily found using Pythogorus:
$Ft = \sqrt{Fr^2 +Fz^2}$
And the component (Fc) of the total force that points directly at the centre is found using:
$Fc = \cos\left( \tan^{-1}\left(\frac{Fz}{Fr}\right) +
\tan^{-1}\left(\frac{-r}{z}\right) - \pi/2 \right) Ft $
$= \cos\left( \tan^{-1}\left(\frac{Fz}{Fr}\right) +
\tan^{-1}\left(\frac{-r}{\sqrt{(a^2 - a^2 e^2) (1 - r^2/a^2)}}\right) - \pi/2 \right) Ft $
Below is a plot of Fr, Fz, Ft, FC and z rom r=0 to r=a:.
This plot is colour coded the same as the force vectors in the diagram above it. The magenta curve (z) is the elliptical profile of the oblate at e = 0.99. The gold curve (Fz) is the vertical component of the surface force and drops off to zero going towards the equator. The blue curve is the horizontal component (Fr) that points at the vertical z axis. This component steadily increases gong from he centre to the equator. The green curve (Ft) is the resultant for Fr and Fz and only drops off by about 20% going from the centre out to the equator. The red curve (Fc) is the component of the total surface force that points directly at the centre and after an initial drop off starts to increase linearly going outwards.
If we naively assumed the central force was proportional to $GMm/R^2$ (where R = $\sqrt{r^2+z^2}$ ), then the required balancing centrifugal force would require a velocity of $u = \sqrt{\frac{G M m}{R}}$, while using the central pointing force Fc derived above, the required velocity to produce the required balancing centrifugal force is given by $v = \sqrt{ Fc \ \ R}$ and these two velocity curve profiles are plotted below for comparison:
The classic assumption (blue) for stable orbital velocities, drops off with increasing distance from the centre, while the exact solution presented here, shows the required balancing velocities increase with increasing distance from the centre. As a side note, it is worth noting that orbits around a galaxy spend 50% of their time above the equatorial plane and the other 50% below the equatorial plane and there is plenty of potential for collisions at the cross over points.
As a sanity check it is possible to find the solution using the surface of equal force (as described in the paper) and assuming the resultant force is normal to this equal force surface and the results are similar for direction.
A typical elliptical galaxy with a central bulge could be modelled by adding the above equations for the flattened disk, to the usual $F = GMm/r^2$ to allow for the extra mass of the bulge at the centre modelled as a sphere. By using extreme values approaching unity for the eccentricity, a model of an almost flat disk can be approximated, which is otherwise very difficult to find an analytical solution for.