For the counterexample, where f is continuous but the second derivatives of [itex]\phi[/itex] blow up. This will occur, e.g., on the edges of a uniform density cube and also if the density in the cube drops off inversely proportional to 1/log(distance to surface of cube), which is continuous. More generally you can replace 'cube' with any solid body whose surface has ridges.
The details are rather involved, but I'll try my best. I don't know a simple and quick method. First we can derive a general expression for the second derivatives. The following derivative is easily calculated for x != 0. I'm using [itex]\delta_{ij}[/itex] for the (nonDirac) delta function equal to 1 if i=j and 0 otherwise.
[tex]
\nabla_i\nabla_j (1/x) = \nabla_i(\hat x_j/x^2) = (3\hat x_i\hat x_j\delta_{ij})/x^3.
[/tex]
Then, the derivatives in the sense of distributions (including point at 0) can be calculated using the divergence theorem as
[tex]
\nabla_i\nabla_j (1/x) = \frac{4}{3}\pi\delta_{ij}\delta(x)+\lim_{r\to 0}1_{\{x>r\}}(3\hat x_i\hat x_j\delta_{ij})/x^3.
[/tex]
This is an equality of distributions, meaning that you integrate vs a smooth function. The limit r>0 on the rhs is similarly a limit in distribution so that the limit is taken after the integration. You can't just set r to 0, because that term wouldn't be locally integrable.
Plugging this into the definition of [itex]\phi[/itex] gives the following
[tex]
\nabla_i\nabla_j\phi(x)=\frac{1}{3}\delta_{ij}f(x)+\lim_{r\to 0}\frac{1}{4\pi}\int_{y>r} f(x+y)\frac{\delta_{ij}3\hat y_i\hat y_j}{y^3}\,dy.
[/tex]
Letting [itex]d\sigma(y)[/itex] be the area integral on S  the sphere of radius 1 centered at 0, we can change variables.
[tex]
\nabla_i\nabla_j\phi(x)=\frac{1}{3}\delta_{ij}f(x)+\lim_{r\to 0}\frac{1}{4\pi}\int_{r}^\infty \int_S f(x+ty)(\delta_{ij}3 y_i y_j)\,d\sigma(y)\,\frac{dt}{t}.
[/tex]
The integral dt/t will blow up as the lower limit r>0. However, the integral of [itex]\delta_{ij}3y_iy_j[/itex] over the sphere S is zero, so if you Taylor expand f about x, the zeroth order term will drop out of the integral above and it remains finite as r>0. Actually, the integral of [itex]\delta_{ij}3y_iy_j[/itex] is also zero over a hemisphere, so the same applies if x is on a smooth boundary surface between two regions where f is smooth. E.g. the derivatives of the gravitational field don't blow up on the boundary of a uniform density ball.
However, it would blow up at the edges of a uniform density cube. Suppose that x is on such an edge with the x_{1} direction pointing inwards and bisecting the angle of the two adjacent faces.
[tex]
\left(\frac{\partial}{\partial x_1}\right)^2\phi(x)=\frac{1}{3}f(x)+\lim_{r\to 0}\frac{1}{4\pi}\int_{r}^\infty \int_S f(x+ty)(13 y_1^2)\,d\sigma(y)\,\frac{dt}{t}.
[/tex]
For small t, f(x+ty) will then be nonzero on a 90 degree 'wedge' on the surface of the sphere [itex]y\in S[/itex], and for y_{1}>= 0. This will cause the integral of f(x+ty)(13y_{1}^{2}) to be strictly negative and nonvanishing as t>0.
So the integral above will be negatively proportional to [itex]\int_{r}^\cdot dt/t[/itex] in the limit as r>0, which blows up at rate log(r).
Finally suppose that the density f inside the cube drops off at rate 1/log(u) towards the surface of the cube where u = distance to edge. Then, f is continuous, and the integral of f(x+ty)(13y_{1}^{2}) will be negative and going to zero at rate 1/log(t) as t> 0. So, the integral above will be negatively proportional to [itex]\int_r^\cdot dt/t\log t[/itex] as r>0, which blows up at rate log(log(r)). The second derivative above will diverge to negative infinity.
