The first step is to find the partial derivative of $$P = \frac{3 M}{(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] \ \ $$ with respect to r, 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.

The result occupies several pages, so I will not post the result here, but it can be seen in this [Mathematica worksheet][1]. 


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.

We can now plot the horizontal gravitational force as a function of r: 

[![enter image description here][2]][2]

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 experimentations shows that the straight line portion of the graph for $F_r$ is approximated fairly well by the following:

$$F_r = \frac{r}{a} \left( \frac{M}{a^2} \exp{(e^2)}^{\left(\frac{2 e^2}{3}\right)} + \frac{1}{\sqrt{1 - e^2}} \right)$$  

The important thing to note is that everything in the equation is a constant except for r and the equation is 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$ 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 =  \sqrt{a^2-r^2} \left(\frac{M}{a^3} (\exp{e^2})^{\left(\frac{9 e^2}{10}\right)}  + \frac{1}{a \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.    

[![enter image description here][3]][3]

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.  



  [1]: https://www.wolframcloud.com/obj/kpegrume/Published/OblateForce.nb
  [2]: https://i.sstatic.net/zOJdk.jpg
  [3]: https://i.sstatic.net/2bWhs.jpg