In principle I don't think your question is really any different than an electrostatics problem where you have some charge distribution in space here and you want to find out what the field strength is over there. The only difference is that we're talking about a mass distribution and not a charge distribution. At first, I was first thinking you could use Gauss' Law, but then I realized that wouldn't work because the field at the surface actually changes in magnitude depending on where you are. I think the hard way to do this is to use Coulomb's Law and just grind through the math.
[itex]\Phi(r,\theta,\phi) = ∫\frac{\rho(\vec{r}')}{|r-r'|} d^3r'[/itex]
The tricky part is doing the actual integral. I think the easiest way might be to parametrize the ellipsoid then do it in spherical coordinates, but that seems like a HUGE pain because I'm not sure how to handle the denominator. The parameters I was trying to use: [itex]\hat{r} = a\sin\theta\cos\phi \hat{x} + b\sin\theta\sin\phi \hat{y} + ccos\theta \hat{z}[/itex] where [itex]\vec{r} = r\hat{r}, \vec{r'} = r\hat{r}'[/itex], and [itex]\rho(\vec{r}') = \rho_0[/itex]
I was also thinking maybe you could use the addition theorem for spherical harmonics, but I'm not sure if that applies here since the parametrization for an ellipse is slightly different than for a sphere. Maybe the integral would be easier to handle in elliptical coordinates?
You could also get an approximate solution for the field by doing a multiple expansion.