Imagine a long bar floating in space. What force does it exert on itself in the middle due to gravity? Problem
If you had a long bar floating in space, what would be the compressive force at the centre of the bar, due to the self-weight of both ends?
Diagram - what is the force at point X in the middle of the bar?:
<----------------------L--------------------->, total mass M
=======================X======================   <- the bar
                 F---> X <---F

Summary
You should be able to simplify by cutting the bar into pieces, but that gives a different answer depending on how many pieces you use (see below). So the simplification must be wrong - but why?
My approach
Split bar in two
So, one approximation would be to cut the bar in half - two pieces of length L/2, mass M/2:
         (M/2)<-------L/2------->(M/2)
           #1          X           #2   <- bar approximated as blobs #1 and #2

Force at X is G(M1.M2)/(R^2) = G (M/2)^2 / (L/2)^2 = G M^2 / L^2
Or Fx / (G. M^2 / L^2) = 1
But is that really valid? If so, shouldn't you get the same answer if you split the bar into four pieces?
Split bar into four
   (M/4)<-L/4->(M/4)<-L/4->(M/4)<-L/4->(M/4)
     #1          #2    X     #3          #4

My assumption is that the force at X is the sum of the attractions of each blob on the left to every blob on the right.
Force at X = #1<>#3 + #1<>#4 + #2<>#3 + #2<>#4
('<>' being force between blobs #x and #y).
Fx / (G.M^2 / L^2) = (2/4)^-2 + (3/4)^-2 + (1/4)^-2 + (2/4)^-2 = 1.61
This is bigger than the previous result (1.61 vs 1).
Split bar into six
Similarly, if you split into 6 blobs, the total force comes out as:
Fx / (G.M^2 / L^2) = (3/6)^-2+(4/6)^-2+(5/6)^-2 + (2/6)^-2+(3/6)^-2+(4/6)^-2 + (1/6)^-2+(2/6)^-2+(3/6)^-2
Fx / (G.M^2 / L^2) = 2.00
So what's wrong with my approach? And what is the real answer?
So it seems the more pieces we split the bar into, the larger the result gets. There's clearly something wrong with my assumptions! - but what?
I'd be very glad if someone here could explain this. Thanks!
EDIT
As Peter Shor pointed out, my calculations had some dodgy algebra and I'd calculated $$L^2/M^2$$ values rather than $$M^2/L^2$$. I've now corrected that - the value still increases as you divide into more masses.
I'll do a bit more work with more divisions and see if this leads to convergence or not.
 A: The reason you encountered higher and higher pressure at the center of the rod as you cut it into more pieces is that you were essentially approximating an integral, but the integral diverges ("is infinity" colloquially).  When the rod has zero thickness, but still has mass, the density of the matter is infinite, and this leads to infinitely strong gravitational forces.
To answer this question, we'll imagine the cylinder has some small, finite radius $R$.  We want to find the force between the two halves of the cylinder.  We'll let one half just sit stationary in space.  It will create a gravitational potential.  Then we'll grab the other half and pull it away to some distance $d$.  The gravitational potential energy is a function of $d$.  The force between the two halves of the cylinder is the derivative of the gravitational potential energy with respect to $d$ when $d=0$.
The problem described above is too hard.  It is quite difficult to calculate the gravitational potential of a cylinder at an arbitrary point.  The gravitational potential of a point mass is just $-Gm/r$, but for a cylinder that extends out in three dimensions, we need to replace $m$ with the density $\rho$ and then integrate over the mass of the entire cylinder.  The expression for $r$, the distance from an arbitrary point outside the cylinder to a point inside it, is not very tractable.  
However, at a point on the axis of the cylinder, the gravitational potential is more accessible due to the extra symmetry.  If we set up cylindrical coordinates with the axis of the cylinder along the z-axis, and then integrate over the bottom half of the cylinder, we get
$V(z) = \int_{z'=0}^{-L/2}\int_{r=0}^{R}\int_{\theta=0}^{2\Pi} \frac{G\rho}{\sqrt{(z-z')^2+r^2}} r\textrm{d}\theta\textrm{d}r\textrm{d}z'$
and doing the integral of $\theta$ it's
$V(z) = 2\pi G\rho\int_{z'=0}^{-L/2}\int_{r=0}^{R} \frac{1}{\sqrt{(z-z')^2+r^2}} r\textrm{d}r\textrm{d}z'$.
This allows us to make an approximation.  Although the half of the cylinder we use to calculate the potential must have finite width, we can calculate the potential energy by assuming that the other half of the cylinder is located perfectly along the axis.  As long as the radius of the cylinder is very small compared to the length, this is a valid approximation.  So the potential energy comes from integrating the previous expression for $V$ along the $z$-axis for the length of the cylinder.
We don't actually want the potential energy, but the derivative of the potential energy.  So we imagine moving the top half up the cylinder up a little bit $dz$, and ask how the potential energy changes.
Moving the entire top half of the cylinder up by $dz$ is equivalent to taking a piece of thickness $dz$ and slicing it off the bottom and moving it to the top.  So we really just need to find the difference in the potential between the top and bottom of the top half of the cylinder and multiply by the mass-per-unit-length of the cylinder.
The force between the two halves of the cylinder is $\frac{M}{L}[V(L/2) - V(0)]$
That still leaves two integrals to evaluate.  $V(L/2)$ is easy, because it's far away from the half of the cylinder providing the gravitational potential (compared to $R$).  That lets us approximate
$V(L/2) = \frac{-GM}{L} \int_{-L/2}^0 \frac{1}{L/2-x}dx = -\frac{GM}{L}\ln 2$.
The integral for $V(0)$ is trickier, so I put it in Mathematica and got
$V(0) = -\frac{GM}{L}\textrm{arcsinh}\left(\frac{L}{2R}\right)$.  
In the regime we are interested in ($R$ small compared to $L$) the $\sinh(x)$ is just $e^x/2$, so this simplifies to
$V(0) = -\frac{GM}{L} \ln\left(\frac{L}{R}\right)$
This gives a final answer for the force
$F = \frac{M^2G}{L^2}\ln\left(\frac{L}{2R}\right)$
A: As other people answered, the force is a vector quantity, so the net force at the center is indeed zero. But it seems what you have in mind is actually the pressure (stress) at that particular point. 
Firstly, regarding your question. What you are computing is
$$ \sum_{\substack{i\in\text{first part}\\j\in\text{second part}}} \mathbf F_{ij} $$
the total force acting on each particle on the first part due to the second part. This is not equivalent to "force at X". But anyway, but why do you get infinity? Because this is the answer.
If we compute the gravitational field due to the 2nd part ($\lambda=M/L$):
$$ g(h) = G\lambda \int_0^{L/2} \frac{1}{(z+h)^2}dz = \frac{G\lambda L}{h(2h+L)}, $$
we find that the gravitational force scales as $\frac1h$. Now suppose you split the first part into $N$ blobs. The position of the $i$-th blob closest to the center will be $h_i=\frac{Li}{2N}$, its mass will be $\frac{\lambda L}{2N}$, so together the contribution of force of each blob will be $F_i = m_ig(h_i) \sim \frac{1/N}{i/N} = \frac1i$, and the total force is $\sim\sum_{i=1}^N\frac1i\sim\ln N$ which diverges logarithmically to $\infty$.
This sounds unphysical, but this is because the setting itself is unphysical — you can't have a mass with no "width" but finite mass, this makes the density infinite, and this causes the infinite force as the two part are right next to each other.
So let's makes sense when the bar has size. So assume that bar is a cylinder of radius $\epsilon \to 0$. Also assume the bar's density is $\rho = M/\pi\epsilon^2L$. The gravity field along the z direction is then
$$ g(h) = G\rho \int_0^{L/2} \int_0^{2\pi} \int_0^{\epsilon} \frac{h+z}{((h+z)^2+r^2)^{3/2}} rdrd\phi dz = G\rho\pi  \left(2 \sqrt{h^2+\epsilon ^2}-\sqrt{(2 h+L)^2+4 \epsilon ^2}+L\right), $$
and the force is now bounded above by a constant $G\rho\pi(2\epsilon+L-\sqrt{4\epsilon^2+L^2})$ for varying $h$, and the total force by partitioning the 1st part would be $\sim \sum_{i=1}^N\frac{\text{constant}}N = \text{constant}$ as you would expect. (Nevertheless, this constant will diverge with $1/\epsilon$ as $\epsilon\to0$.)

But the result above has nothing to do with the pressure. At precisely the center, there is no force, and the magnitude of force scales as $r$ in its immediate neighborhood. 
Let's mark a small square membrane of size $a\times a\times dz$ with negligible mass $a^2\rho dz$ at position $z\to0$ above the center of the bar to record the pressure.
The membrane will experience a net push from above. Since the membrane is small, we could assume it is a uniform force. The total force that pushes the membrane downwards is
$$\begin{aligned} F &= a^2\rho dz g_{\text{total}}(z) \\
&= a^2\rho dz\cdot G\rho\pi \left(\sqrt{(L+2z)^2+4\epsilon^2}-\sqrt{(L-2z)^2+4\epsilon^2}\right) \\
&\sim  \frac{4a^2z\rho dz\cdot G\rho L\pi }{\sqrt{L^2+4\epsilon^2}}
 \end{aligned}$$
and thus the pressure is 
$$ dP = -\frac F{a^2} = -\frac{4zG\rho^2 L\pi dz }{\sqrt{L^2+4\epsilon^2}} \sim 4zdz \frac{GM^2}{L^2\epsilon^4} $$
$$ \implies P \sim -z^2 \frac{GM^2}{2L^2\epsilon^4} $$
note that the $z$ is still here. We see that the pressure is finite and maximal when we measure it precisely at the center (and decreases away from it). But this will diverge if we let $\epsilon\to0$.
A: For one thing, you've made a careless algebra mistake.
You're computing $(L/M)^2$ instead of $(M/L)^2$.
