Determining geodesics between two points in curved space In order to determine the geodesics between two points, one must solve the geodesic differential equations, which are as following
\begin{align}
p &= u'(s)\\
q &= v'(s)\\
p' + \Gamma^0_{00}p^2 + 2\Gamma^0_{01}pq + \Gamma^0_{11}q^2 &= 0\\ 
q' + \Gamma^1_{00}p^2 + 2\Gamma^1_{01}pq + \Gamma^1_{11}q^2 &= 0
\end{align}
where $\Gamma^i_{jk}$ is the Christoffel symbol of second kind, and the initial conditions are $u(s_0) = u_0,\, \, v(s_0) = v_0,\, \,  u'(s_0) = p_0,\, \, v'(s_0) = q_0$.
In order to solve these equations, one can only provide one point $(u,v)$, and the direction $(u',v')$. How is it possible to determine the geodesics between two points ?
For a sphere, both points will lie on the great circle. So here it is sufficient to provide the initial point and direction. But is it possible to accomplish this by only providing both points, if information such as direction is not known a priori ?
 A: This is not always possible, consider for example the punctured (at $(0,0)$) plane with the usual flat metric, and the points $(-1, 0)$ and $(1, 0).$
A: I assume that   imranal     asks  how to  find numerically a geodesic connecting two given points if the  connection is given. 
One way to do it  is to implement the solution of the ODE system he wrote in his question numerically (there are many effective ways for it) and then use this implimintation  to find the correct initial velocity vector $(u', v')$  such that the geodesic starting from $(u, v)$ in the direction $(u', v')$ hits the second point.  One can  slighly improve this search but the principle remains the same. 
Alternative method  which actually needs that the connection you have comes from a Riemannian metric 
would be to use the curve shortening flow: start with any curve connecting these two points (say, a straight line)  and then
 start to deform  the curve  such that the velocity vector of the deformation is orthogonal to the velocity vector of the curve and its length is    the  geodesic curvature times some fixed function which vanishes at the first and  last points of the curve. 
This deformation  should  normally converge to a geodesic. 
I expect that there should be  effective implimentation of this  idea as a numerical scheme but I am not an expert in these questions. 
Of course geodesic does not always exist  (as explained  by Igor) and   may be is not unique (as explained by Narasimham) or/and the second procedure may converge to a geodesic which is not a minimal one. 
A: If two points $P,Q$ which are not umbilical are given, then there are in general several geodesics but only one that is shortest among them all. Their  $ u^{'}, v^{'}$ are different, but unique for any one of the geodesic choices.
If $P,u^{'}, v^{'}$ are given there is no guarantee that the line so defined passes through $Q.$ 
These can be physically verified between two points on a cylinder/cone using a taut thread for example. 
EDIT1:
It can be compared to the dynamic geodesic trajectory situation when a gun is fired with a given velocity and inclination from start point $P$. The initial velocity and direction should be adjusted to make trajectory pass through another desired point $Q$ in the vertical plane, as it cannot meet a desired point with arbitrary setting of dynamic parameters like $ u^{'}, v^{'} $ in case of geometrical geodesics.
EDIT2:
The solution also reflects in numerical procedures adopted. We have both initial value and boundary value problems. The former is straight forward and for the latter position/ derivative values at $P$ and $Q$ are given as input and numerical iterative (again!) shoot-through technique is employed.
EDIT3:
Actually the answer to your question is determination of the common geodesic invariant through  $P,Q$. Geodesic curvature $k_g$  vanishes for all geodesics, so condition for them to be on the same geodesic line is that they be on same filament of the tangent fibre bundle  sharing the same invariant.
The parameters of the two points go to determine it..
For a surface of revolution the condition leads to Clairaut's Law.. a particularly easy relation. Geodesic constant condition  is derivable as Liouville's relation generally for the $k_g=0$ condition from Christoffel symbols.
