On $\mathbb{R}^3$ with coordinates $x,y,z$ consider the Riemannian metric $g=\displaystyle{\frac{dx^2+dy^2+dz^2}{x^2+y^2}}$ defined on $X:=\mathbb{R}^3\setminus\{(0,0,z)\}$.
Given any point $p_1\in X$ and $p_2$ any point in a sufficiently small neighborhood of $p_1$, I want to compute the distance $d(p_1,p_2)$ with respect to the metric $d$ which is the intrinsic metric induced by $g$.
Take for example $p_1=(1,1,1)$ and $p_2=(3/2,-3/2,2)$ how do I compute $d(p_1,p_2)?$
The only way I've figured to compute $d(p_1,p_2)$ is to compute the geodesic $\gamma$ from $p_1$ to $p_2$ and then to compute its length. I've computed the geodesic equations of $g$ ($\gamma_1$ is the component of $\gamma$ in the coordinate $x$ and so on):
$$\ddot\gamma_1=\displaystyle{\frac{2y\dot\gamma_1\dot\gamma_2+x(\dot\gamma_1^2-\dot\gamma_2^2-\dot\gamma_3^2)}{x^2+y^2}}$$
$$\ddot\gamma_2=\displaystyle{\frac{2x\dot\gamma_1\dot\gamma_2-y(\dot\gamma_1^2-\dot\gamma_2^2+\dot\gamma_3^2)}{x^2+y^2}}$$
$$\ddot\gamma_3=\displaystyle{\frac{2(x\dot\gamma_1+y\dot\gamma_2)\dot\gamma_3}{x^2+y^2}}$$
But I don't know how to integrate these differential equations to get the equation of the geodesic from $p_1$ to $p_2$ and then compute its length.
Note that your metric depends only on $x^2 + y^2$ so it makes sense to exploit this symmetry by performing a change of coordinates to cylindrical coordinates:
$$ x = \rho \cos \theta, y = \rho \sin \theta, z = z. $$
In those coordinates, the metric has the simpler form
$$ g = \frac{1}{\rho^2} (d\rho^2 + dz^2) + d\theta^2. $$
Working with the order $(\rho, z, \theta)$ for the coordinates and assuming I haven't done a computation mistake, the geodesic equations for $\gamma = (\gamma_1,\gamma_2,\gamma_3)$ (where $\gamma_1$ corresponds to $\rho$, $\gamma_2$ corresponds to $z$, etc) take the form
$$ \ddot{\gamma}_1 - \frac{1}{\gamma_1} \left( \dot{\gamma}_1^2 - \dot{\gamma}_2^2 \right) = 0, \\ \ddot{\gamma}_2 - \frac{2}{\gamma_1} \left( \dot{\gamma}_1 \dot{\gamma}_2 \right) = 0, \\ \ddot{\gamma}_3 = 0. $$
The equation for $\gamma_3$ can be solved directly and the second equation reveals the relation $\dot{\gamma_2} = A \gamma_1^2$ which can be used to solve completely the geodesic equations.