I am not sure what is the path $C$ you are integrating over? In your definition you evaluate $U(C)$ which in the present case of force is independent on the explicit path you choose but still depends on initial and final point, i.e. $U(p_1,p_2)$. In your final result it seems you are actually 'walking' three times the path $p_1=(-\infty,y,z)$ to $p_2=(-x,y,z)$. only by relabeling your coordinate system. Thus the factor of $3$ in your final formula.

Best practice is to parametrize the path $C$ in terms of a function $\vec{s}(\lambda)=(x(\lambda),y(\lambda),z(\lambda))$ and integrating over $\lambda$.

In Detail

In order to evaluate the potential you have to know the path you are integrating over: $$ U(C)=\int_{\vec{r}\in C} d\vec{r}\cdot \vec{F}(\vec{r})$$ In the present case the function (the force) is conservative (vanishing curl) and thus the result only depends on the initial and final point $\vec{r}_{i}$ and $\vec{r}_f$ of $C$. Now, what is the initial point and the final point in your setting. You say $C$ should connect infinity with the point $\vec{r}_f=(x,y,z)$. Thus, what is $\vec{r}_i$ here? Basically you can choose any 'boundary' point of $\mathbb{R}^3$. A good choice would be $\vec{r}_i=(-\infty,y,z)$, though other choices as for instance $\vec{r}_i=(-\infty,-\infty,-\infty)$ are also suitable. Now you can basically choose any path connecting $\vec{r}_i$ and $\vec{r}_f$ it will give the same result. So take a very simple one $C_s$: $\vec{s}(\lambda)=(\lambda ,y,z)$ with $\lambda\in (-\infty,x]$. The integral is then parametrized and what you have to evaluate is

$$ U(\vec{r}_i,\vec{r}_f)=\int_{\vec{s}\in C_s} d\vec{r}\cdot \vec{F}(\vec{r})\equiv \int_{-\infty}^{x}d \lambda \frac{d\vec{s}(\lambda)}{d\lambda}\cdot \vec{F}(\vec{s}(\lambda))=\int_{-\infty}^{x}d \lambda \, \hat{e}_x\cdot \vec{F}(\vec{s}(\lambda))$$ Now insert the path $\vec{s}(\lambda)=(\lambda ,y,z)$ and you obtain $$U(\vec{r}_i,\vec{r}_f)= \int_{-\infty}^{x}d \lambda \, \left(-G m M \frac{\lambda}{(\lambda^2+y^2+z^2)^{3/2}} \right) $$ Other paths will yield the same result, but in your case you follow a route starting from three different points $(-\infty,y,z)$, $(x,-\infty,z)$, and $(x,y,-\infty)$ and terminate at $(x,y,z)$. Thus factor of $3$.

Surely, you can convert the result into spherical coorindates and use a path that is radial.

I am not sure what is the path $C$ you are integrating over? In your definition you evaluate $U(C)$ which in the present case of force is independent on the explicit path you choose but still depends on initial and final point, i.e. $U(p_1,p_2)$. In your final result it seems you are actually 'walking' three times the path $p_1=(-\infty,y,z)$ to $p_2=(-x,y,z)$. only by relabeling your coordinate system. Thus the factor of $3$ in your final formula.

Best practice is to parametrize the path $C$ in terms of a function $\vec{s}(\lambda)=(x(\lambda),y(\lambda),z(\lambda))$ and integrating over $\lambda$.

In Detail

In order to evaluate the potential you have to know the path you are integrating over: $$ U(C)=-\int_{\vec{r}\in C} d\vec{r}\cdot \vec{F}(\vec{r})$$ In the present case the function (the force) is conservative (vanishing curl) and thus the result only depends on the initial and final point $\vec{r}_{i}$ and $\vec{r}_f$ of $C$. Now, what is the initial point and the final point in your setting. You say $C$ should connect infinity with the point $\vec{r}_f=(x,y,z)$. Thus, what is $\vec{r}_i$ here? Basically you can choose any 'boundary' point of $\mathbb{R}^3$. A good choice would be $\vec{r}_i=(-\infty,y,z)$, though other choices as for instance $\vec{r}_i=(-\infty,-\infty,-\infty)$ are also suitable. Now you can basically choose any path connecting $\vec{r}_i$ and $\vec{r}_f$ it will give the same result. So take a very simple one $C_s$: $\vec{s}(\lambda)=(\lambda ,y,z)$ with $\lambda\in (-\infty,x]$. The integral is then parametrized and what you have to evaluate is

$$ U(\vec{r}_i,\vec{r}_f)=-\int_{\vec{s}\in C_s} d\vec{r}\cdot \vec{F}(\vec{r})\equiv -\int_{-\infty}^{x}d \lambda \frac{d\vec{s}(\lambda)}{d\lambda}\cdot \vec{F}(\vec{s}(\lambda))=-\int_{-\infty}^{x}d \lambda \, \hat{e}_x\cdot \vec{F}(\vec{s}(\lambda))$$ Now insert the path $\vec{s}(\lambda)=(\lambda ,y,z)$ and you obtain $$U(\vec{r}_i,\vec{r}_f)= -\int_{-\infty}^{x}d \lambda \, \left(-G m M \frac{\lambda}{(\lambda^2+y^2+z^2)^{3/2}} \right) $$ Other paths will yield the same result, but in your case you follow a route starting from three different points $(-\infty,y,z)$, $(x,-\infty,z)$, and $(x,y,-\infty)$ and terminate at $(x,y,z)$. Thus a factor of $3$ appears.

The same holds true for the initial point $\vec{r}_i=(-\infty,-\infty,-\infty)$. Basically, you can think of any path connecting $\vec{r}_i$ with $\vec{r}_f$ and then try to find a function that describes this curve. A naive and very simple choice would be a straight line: $ \vec{s}(\lambda)=(\lambda x, \lambda y, \lambda z)$ with $\lambda\in (-\infty,1]$ and $\frac{d\vec{s}(\lambda)}{d\lambda}=x\hat{e}_x+y\hat{e}_y+z \hat{e}_z$. However, actually you are crossing the origin for which the force diverges.

Nevertheless, just for simplicity we following this path along which the force assumes the form $$\vec{F}(\vec{s}(\lambda))=-GmM \frac{ x\lambda \hat{e}_x+y\lambda \hat{e}_y+z\lambda \hat{e}_z}{( (x\lambda)^2+(y\lambda)^2+(z\lambda)^2 )^{3/2}}= -GmM \frac{ x\hat{e}_x+y\hat{e}_y+z \hat{e}_z}{( x^2+y^2+z^2 )^{3/2}} \frac{1}{\lambda^2}$$ Hence, the kernel of the integral looks like: $$\frac{d\vec{s}(\lambda)}{d\lambda}\cdot \vec{F}(\vec{s}(\lambda))= -GmM \frac{x^2+y^2+z^2}{(x^2+y^2+z^2)^{3/2}}\frac{1}{\lambda^2}\equiv -\frac{GmM}{r} \frac{1}{\lambda^2} $$ whereby $r=\sqrt{x^2+y^2+z^2}$ the radial of the end-point.

Now, the integral can easily be evaluated: $$U(\vec{r}_i,\vec{r}_f)= -\int_{-\infty}^{1}d \lambda \, \left(- \frac{G mM}{r} \frac{1}{\lambda^2} \right) =-\frac{G mM}{r} \left[\frac{1}{\lambda} \right]^{1}_{-\infty}=-\frac{G mM}{r}$$

Surely, you can convert the result into spherical coordinates and use a path that is radial.

I am not sure what is the path $C$ you are integrating over? In your definition you evaluate $U(C)$ which in the present case of force is independent on the explicit path you choose but still depends on initial and final point, i.e. $U(p_1,p_2)$. In your final result it seems you are actually 'walking' three times the path $p_1=(-\infty,y,z)$ to $p_2=(-x,y,z)$. only by relabeling your coordinate system. Thus the factor of $3$ in your final formula.

Best practice is to parametrize the path $C$ in terms of a function $\vec{s}(\lambda)=(x(\lambda),y(\lambda),z(\lambda))$ and integrating over $\lambda$.

I am not sure what is the path $C$ you are integrating over? In your definition you evaluate $U(C)$ which in the present case of force is independent on the explicit path you choose but still depends on initial and final point, i.e. $U(p_1,p_2)$. In your final result it seems you are actually 'walking' three times the path $p_1=(-\infty,y,z)$ to $p_2=(-x,y,z)$. only by relabeling your coordinate system. Thus the factor of $3$ in your final formula.

Best practice is to parametrize the path $C$ in terms of a function $\vec{s}(\lambda)=(x(\lambda),y(\lambda),z(\lambda))$ and integrating over $\lambda$.

In Detail

In order to evaluate the potential you have to know the path you are integrating over: $$ U(C)=\int_{\vec{r}\in C} d\vec{r}\cdot \vec{F}(\vec{r})$$ In the present case the function (the force) is conservative (vanishing curl) and thus the result only depends on the initial and final point $\vec{r}_{i}$ and $\vec{r}_f$ of $C$. Now, what is the initial point and the final point in your setting. You say $C$ should connect infinity with the point $\vec{r}_f=(x,y,z)$. Thus, what is $\vec{r}_i$ here? Basically you can choose any 'boundary' point of $\mathbb{R}^3$. A good choice would be $\vec{r}_i=(-\infty,y,z)$, though other choices as for instance $\vec{r}_i=(-\infty,-\infty,-\infty)$ are also suitable. Now you can basically choose any path connecting $\vec{r}_i$ and $\vec{r}_f$ it will give the same result. So take a very simple one $C_s$: $\vec{s}(\lambda)=(\lambda ,y,z)$ with $\lambda\in (-\infty,x]$. The integral is then parametrized and what you have to evaluate is

$$ U(\vec{r}_i,\vec{r}_f)=\int_{\vec{s}\in C_s} d\vec{r}\cdot \vec{F}(\vec{r})\equiv \int_{-\infty}^{x}d \lambda \frac{d\vec{s}(\lambda)}{d\lambda}\cdot \vec{F}(\vec{s}(\lambda))=\int_{-\infty}^{x}d \lambda \, \hat{e}_x\cdot \vec{F}(\vec{s}(\lambda))$$ Now insert the path $\vec{s}(\lambda)=(\lambda ,y,z)$ and you obtain $$U(\vec{r}_i,\vec{r}_f)= \int_{-\infty}^{x}d \lambda \, \left(-G m M \frac{\lambda}{(\lambda^2+y^2+z^2)^{3/2}} \right) $$ Other paths will yield the same result, but in your case you follow a route starting from three different points $(-\infty,y,z)$, $(x,-\infty,z)$, and $(x,y,-\infty)$ and terminate at $(x,y,z)$. Thus factor of $3$.

Surely, you can convert the result into spherical coorindates and use a path that is radial.