Finding The Maximal Gravity Here is my question

Given a point $P$ in space, and given a piece of malleable material of
  constant density, how should you shape and place the material in order
  to create the largest possible gravitational field at $P$?

Label the points on the surface by their distance $r$ from $P$, and by the angle $θ$
that the line of this distance subtends with the $x$-axis. Then a small mass $dm$ on
the surface provides an $x$-component of the gravitational field equal to
$$F_{x} = \cos θ \frac {G\; dm}{r^2} $$
Any help is appreciated!
 A: Before we start, let us rephrase the question to a way I understand:

Given a point $P \in \mathbb{R}^3$ and a piece of malleable material of mass $M$, constant density $\rho$, determine the shape of the material (i.e. a region $\Omega \subset \mathbb{R}^3$ with volume $V = \frac{M}{\rho}$) so that the 
  gravitation force it asserted to a test charge of mass $m$ at $P$
  $$\vec{F} = Gm\rho \int_{\Omega} \frac{\vec{x} - \vec{P}}{|\vec{x}-\vec{P}|^3} dxdydz$$
  is maximized in magnitude.

The problem has a lot of symmetry, it is invariant under translation of $P$ and rotation of directional axis. WOLOG, we only need to consider the case $P$ is the origin and the coordinate axis is chosen so that the "largest" gravitation force is pointing towards the $z$ direction. i.e.
$$\text{maximize}\;\frac{F_z}{Gm\rho} = \int_{\Omega} \frac{z}{(x^2+y^2+z^2)^{3/2}} dx dy dz
\quad\text{ subject to }\quad \verb/Volume/(\Omega) = V$$
Let $1_\Omega(x,y,z)$ be the indicator function for $\Omega$, it will be a solution to following problem
$$\text{maximize}\;\int f(x,y,z)\frac{z}{(x^2+y^2+z^2)^{3/2}}dxdydz
\quad\text{ subject to }\quad
\begin{cases}\int f(x,y,z) dxdydz = V\\
f(x,y,z) = 0 \text{ or } 1
\end{cases}
$$
By a variant of Rearrangement inequality for integral, the LHS is maximized
when $f$ has the form
$$f(x,y,z) = \begin{cases}
1, & \frac{z}{r^3} > \frac{1}{a^2}\\
0, & \frac{z}{r^3} < \frac{1}{a^2}
\end{cases}
$$
for some suitable chosen $a$. In terms of polar coordinate
$$[0,\infty) \times [0,\pi] \times [ 0,2\pi ] \ni (r,\theta,\phi) \quad\mapsto\quad
(x,y,z) = (r\sin\theta\cos\phi,r\sin\theta\sin\phi,r\cos\theta) \in \mathbb{R}^3
$$
$\Omega$ is the region such that
$$\frac{z}{(x^2+y^2+z^2)^{3/2}} = \frac{\cos\theta}{r^2} \ge \frac{1}{a^2}\quad\iff\quad r \le a\sqrt{\cos\theta}, \theta \in [0,\frac{\pi}{2}]$$
It is a solid of revolution which can be generated by rotating the dipole curve $(x^2+z^2)^3 = z^2$ in $xz$-plane with respect to $z$-axis and then keep the component bounded by the resulting surface above the $xy$-plane.
For any $a$, the volume of $\Omega$ equals to
$$\begin{align}
V = \int_{\Omega} dxdydz 
&= \int_0^{2\pi}\int_0^{\pi/2}\int_0^{a\sqrt{\cos\theta}} r^2 \sin\theta dr d\theta d\phi\\
&= \frac{2\pi a^3}{3} \int_0^{\pi/2}(\cos\theta)^{3/2} \sin\theta d\theta
 = \frac{2\pi a^3}{3} \int_0^1 u^{3/2} du = \frac{4\pi a^3}{15}
\end{align}
$$
This allow us to fix $a$ as $\displaystyle\left(\frac{15 V}{4\pi}\right)^{1/3} = \left(\frac{15M}{4\pi\rho}\right)^{1/3}$. The maximum value of $|\vec{F}|$ satisfies
$$\begin{align}
\frac{|\vec{F}|_{max}}{Gm\rho} = \int_{\Omega} \frac{\cos\theta}{r^2} dxdydz
&= \int_0^{2\pi}\int_0^{\pi/2}\int_0^{a\sqrt{\cos\theta}} \cos\theta \sin\theta dr d\theta d\phi\\
&= 2\pi a \int_0^{\pi/2}(\cos\theta)^{3/2}\sin\theta d\theta
 = 2\pi a \int_0^1 u^{3/2} du = \frac{4\pi a}{5}\\
&= \frac{3V}{a^2}
\end{align}
$$
This leads to $\displaystyle\;|\vec{F}|_{max} = \frac{3GMm}{a^2} = 3GMm\left(\frac{4\pi\rho}{15M}\right)^{2/3}$.
A: Because of axial symmetry, it is enough to determine y(x), 0 ≤ x ≤ D, where diameter D is the farthest point of the surface. Optimal shape means all points of the surface contribute equally to the gravitational pool, so an arbitrary point (x,y) contributes the same as (D,0), therefore y/(x^2+y^2)^(3/2) = 1/D^2. 
Answer: y = x·√(D/x)ⁿ - 1), where n=4/3. (For comparison, equation of a sphere is the same, but with n=1.) My calculations show that compared to a sphere of equal volume, this shape would have a shorter diameter, i.e. if point (0,0) is the "north pole", then distance between two poles is 0.855 of the diameter of a sphere. However, the radius at the equator will be 1.2408 that of a sphere, and the equator itself will be closer to the north pole with thickest point at x=0.4387*D. 
My calculations also show that gravitational pull will be 28% greater than for a sphere of equal volume. 
