Find geodesics given two points

Suppose I want to find the geodesic on a paraboloid going through two points located on the paraboloid, using Mathematica 12.

I chose the paraboloid parametrization as follows:

paraboloid[{u_, v_}] := {Sqrt[u] Sin[v], Sqrt[u] Cos[v], u}

(Source: https://mathworld.wolfram.com/ParaboloidGeodesic.html)

Load the required package to do variational calculus:

Needs["VariationalMethods`"]

Then I set up the Euler-Lagrange equations like so:

eq = EulerEquations[Sqrt[Total[D[paraboloid[{u, v[u]}], u]^2]], v[u],u]

The output looks quite involved, so I was surprised that DSolve could obtain an exact solution using

DSolve[eq, v[u], u]

The two solutions contains ArcTanh, ArcSinh and quite a few roots (one solution below):

v[u_] := (u (ArcTanh[(2 Sqrt[1 + 4 u])/Sqrt[4 - u C[1]]] Sqrt[C[1]] Sqrt[4 - u C[1]] - 
     4 ArcSinh[(Sqrt[1 + 4 u] Sqrt[C[1]])/Sqrt[-16 - C[1]]] Sqrt[-16 -
        C[1]] Sqrt[(4 - u C[1])/(16 + C[1])]))/(Sqrt[C[1]] Sqrt[u^2 (-4 + u C[1])]) + C[2]

Since I am interested in a specific geodesic passing through two given points, e.g. $p_1=(u_1,v_1)=(1,1)$ and $p_2=(u_2,v_2)=(-1,-1)$ on the paraboloid,
I tried to determine the values of C[1] and C[2] like so:

Solve[{v[1] == 1, v[-1] == -1}, {C[1], C[2]}]

Neither Solve nor NSolve would return for at least an hour which is when I stopped the kernel, so how would I go about finding the specific solution making the geodesic go through $(u_1,v_1)$ and $(u_2,v_2)$?

EDIT:
To clarify, I would like to know whether solving for C[1] and C[2] the way I tried is the proper way to obtain the geodesic connecting the points on the given surface. If the principle is sound, how can I let MMA arrive at a (numeric) solution?

UPDATE:
I realized I made a mistake: The paraboloid parametrization at the top makes use of cylinder coordinates, therefore $u$ and $v$ represent the radius and the azimuth angle, respectively. So in order to trace out a complete, non-degenerate paraboloid, we need $u > 0$ and $v in [0,2pi]$. That means the second point through which the specific geodesic should pass cannot, by definition, be $p_2=(-1,-1)$ as above. $(2,frac{3pi}{4})$ is inside the intervals for $u$ and $v$, so that would be OK.

Several posts suggested using FindRoot to solve (systems of) transcendental equations like this one, so I took the solution to the E-L equation from above and modified the parameter list to include the integration constants $C_1$ and $C_2$:

v[u_, C1_, C2_] := (u (ArcTanh[(2 Sqrt[1 + 4 u])/Sqrt[4 - uC1]]Sqrt[C1]Sqrt[4 - u C1] - 4 ArcSinh[(Sqrt[1 + 4 u] Sqrt[C1])/Sqrt[-16 - C1]] Sqrt[-16 - C1] Sqrt[(4 - u C1)/(16 + C1)]))/(Sqrt[C1] Sqrt[u^2 (-4 + u C1)]) + C2

Find $C_1$ and $C_2$ to get the geodesic through $p_1$ and $p_2$:
Let $p_1=(1,0)$, $p_2=(2,frac{3pi}{4})$, so

FindRoot[{v[1, C1, C2] == 0, v[2, C1, C2] == (3 [Pi])/4}, {{C1, 1}, {C2, 1}}]

FindRoot returns {C1 -> -21.7912 - 1.99114*10^-14 I, C2 -> -1.52118 - 0.824159 I}. I suppose I may safely assume the imaginary part of $C_1$ to be zero, but not the one in $C_2$. Also, it says

FindRoot: The line search decreased the step size to within tolerance specified by AccuracyGoal and PrecisionGoal but was unable to ...

I looked at quite a number of network posts that reported the same issue, most of them seem to have specific workarounds which don’t attack the problem simply by increasing WorkingPrecision – how can I get a correct solution?