The problem is that "point masses" are fictitions that make our lives often simpler but only as long as we don't look at the place where the "point mass" sits.
To solve your problem rather consider a mathematically better defined model for the source of your gravitational field first. The most simple you can think of is a sphere with a finite radius ##a## located at the origin of the coordinate system with a constant mass density ##\rho##, i.e.,
$$\rho(\vec{x})=\begin{cases} \rho_0 & \text{for} \quad |\vec{x}| \leq a \\0 & \text{for} \quad |\vec{x}|>a. \end{cases}.$$
Now you can apply the usual math of vector calculus. You know that the gravitational field ##\vec{g}(\vec{x})## is a conservative field, i.e., that there is a scalar potential ##\Phi## such that
$$\vec{g}=-\vec{\nabla} \phi$$
and then there's "Gauss's Law",
$$\vec{\nabla} \cdot \vec{g}=-\Delta \Phi=-4 \pi G \rho,$$
where ##G## is Newton's gravitational constant (the factor ##4 \pi## is conventional).
Now you can solve the problem very easily by first thinking about the symmetries of the problem. Obviously with our choice of a spherically symmetric mass distribution it makes sense to introduce spherical coordinates. Then ##\rho(\vec{x})=\rho(r)## and obviously the ansatz ##\Phi(\vec{x})=\Phi(r)## can be expected to work. Indeed using the known formula for the Laplace operator in spherical coordinates leads to
$$\Delta \Phi(r)=\frac{1}{r} \partial_r^2 (r \Phi)=4 \pi \rho.$$
That's easy to solve for our problem. Let's start with the region ##r>a##. Then the right-hand side vanishes, and we have
$$(r \Phi)''=0 \; \Rightarrow \; r \Phi=A +B r \; \Rightarrow \; \Phi(r)=\frac{A}{r} + B.$$
The constant ##B## is irrelevant and can be set to 0, because adding a constant to the potential doesn't change the field ##\vec{g}=-\vec{\nabla} \Phi##. We can't determine ##A## yet. For that we need to also solve the equation for ##r<a##, where
$$\frac{1}{r} (r \Phi)''=4 \pi \rho_0.$$
Now
$$(r \Phi)''=-4 \pi G \rho_0 r \; \Rightarrow \; (r \Phi)'=2 \pi G \rho_0 r^2 + B' \; \Rightarrow \; r \Phi=\frac{2}{3} \pi \rho_0 r^3 + B' r +A' \; \Rightarrow \; \Phi=-\frac{2}{3} G \pi \rho_0 r^2 + B'+\frac{A'}{r}.$$
Now since there are no singularities in ##\rho## at ##r=0##, there shouldn't be singularities there for ##\Phi## either. Thus we have to set ##A'=0##. Further the potential as well as the field
$$\vec{g}=-\vec{\nabla} \Phi=-\Phi'(r) \vec{e}_r$$
should be continuous at ##r=a##, from which
$$\frac{2}{3} \pi G \rho_0 a^2 +B'=\frac{A}{a}, \qquad (*)\\
\frac{4}{3} \pi G \rho_0 a =-\frac{A}{a^2}.$$
From the latter equation we find
$$A=-\frac{4}{3} \pi G a^3 \rho_0=-M G,$$
where ##M=\rho_0 V=4 \pi a^3 \rho/3## is the total mass of our spherical body.
Finally we get from (*)
$$\frac{M G}{2 a}+B'=-\frac{M G}{a} \; \Rightarrow \; B'=-\frac{3 M G}{2a}.$$
Now let's see what happens, when we make ##a \rightarrow 0## but keeping ##M=\text{const}##. Since ##\rho_0=M/V=3 M/(4 \pi a^3)## we have ##\rho_0 \rightarrow \infty##.
This shows, what's going on for the "point-particle limit". You get
$$\rho(\vec{x})=\begin{cases} 0 & \text{for} \quad \vec{x} \neq 0, \\ \rightarrow \infty & \text{for} \quad \vec{x} \rightarrow 0,
\end{cases}
$$
but
$$M=\int_{V} \mathrm{d}^3 x \rho(\vec{x})$$
for any volume ##V## which contains the origin. Thus the correct limit for the mass density needs an extension of our application of vector calculus to socalled "generalized functions" or "distributions". In this case we get the famous Dirac ##\delta## distribution,
$$\rho(\vec{x})=M \delta^{(3)}(\vec{x}),$$
and
$$\Phi(\vec{x})=-\frac{M G}{r},$$
which diverges for ##r \rightarrow 0##. So you cannot apply the usual notions of vector calculus, where curl, div, grad are given as derivatives with respect to the coordinates for ##r=0##, but everything gets consistent by defining
$$\Delta \Phi(\vec{x})=4 \pi G \rho(\vec{x})=4 \pi G M \delta^{(3)}(\vec{x}),$$
from which you get the extended meaning of the Laplace operator in the sense of generalized functions:
$$\Delta \frac{1}{r} = -4 \pi \delta^{(3)}(\vec{x}).$$
As you see, the apparently "simple" case of a "point mass" needs the introduction of an entire new kind of mathematics, and indeed the Dirac ##\delta## distribution, introduced in a handwaving way by Dirac to make calculations in quantum theory simpler, triggered the development of the entire new mathematical discipline of functional analysis.
