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 M is the mass, "a" is the semi-major axis of the elliptical cross section of the oblate. "e" is the eccentricity of the oblate. "z" is the vertical distance from the equatorial plane and "r" is the horizontal distance from the vertical rotation axis of the point on the surface of the oblate. r and z are cylindrical coordinates but I am ignoring the longitudinal $\phi$ coordinate due the symmetry and r and z are effectively cartesian coordinates on 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.
The exact result for the horizontal component $F_r$ of the surface gravitational force :
$\frac{3 GMm}{ 16 (a e r)^3} \left(4 r^4 (W - 2 Q ) -\frac{2 U r^2 (2 a^2 + r^2 (2 e^2-3) ) }{V} - \frac{W (-a^4 + a^2 r^2 (1 + V) + r^4 S (e^2 - V)) \csc(Q)^2}{V} + \frac{16 r^2 S(a^2 - r^2) \sin(Q)^2}{V W} + \frac{2 r^4 W \sin^2(Q)}{V} \right)$
where $S = (1-e^2)$,
$V = (a^2 - e^2 r^2)/r^2,$
$W = \sqrt{ 2 a^2/r^2 (1 + V) -2 a^4/r^4 + 2 S(e^2 - V)}$
and $Q = 1/2 \arccos(S - a^2/r^2 + V) $.
See this Mathematica worksheet for more details and the calculations.
We can now plot the horizontal gravitational force as a function of r:
In the above plot it can be seen that although the curve of $F_r$ is not trivial, it is in fact a straight line within the confines of $r<a$, i.e. within the radius of the oblate sphere.
A little experimentation shows that the straight line portion of the graph for $F_r$ is approximated fairly well (up to about e=0.999) by the following:
$$F_r = \frac{r G M m}{a} \left( \frac{\exp(2 \ e^4 /3)}{a^2} + \frac{1}{1000 \sqrt{1 - e^2}} \right)$$
where k is just a unit constant in suitable units such as Newtons. The important thing to note is that everything in the equation is a constant except for r and the equation is approximately proportional to $r M$ and not the expected $M/r^2$. It means the surface gravity vector of an oblate sphere that points towards the central axis increases as r gets larger!
Of course the above equation is an approximation, but it is not too difficult to show that for any value of a for the radius of the galaxy and any value of $e$ for the eccentricity of the galaxy, the $r M$ (plus a constant) relationship always holds. For the picky, the exact equations are contained in the Mathematica document linked above.
For completeness, the force component $F_z$ that acts down towards the disk parallel to the rotation z axis, can be approximated by:
$$F_z = \frac{G M m \sqrt{a^2-r^2}}{a} \left(\frac{\exp(9 \ e^4 / 10)}{a^2} + \frac{1}{ 1000 \sqrt{1 - e^2}} \right)$$
With these two approximations we can readily plot the force vectors and visualise the resultant total force vector at the surface and demonstrate that generally, the force vector does not point directly at the centre of the oblate mass.
As a sanity check it is possible to find the solution using the surface of equal force 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 of the disk, a model of an almost flat disk can be approximated, which is otherwise very difficult to find an analytical solution for.