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)$.
This gives a final answer for the force
$F = \frac{M^2G}{L^2}[\textrm{arcsinh}\left(\frac{L}{2R}\right) - \ln 2]$
This is only good in the limit where $R$ is much smaller than $L$. (For instance it gives $F = 0$ for a certain value of $L/2R$, but the approximation is not valid in this regime.)