Calculating time to arbitrary points of distance with initial velocity and non-uniform acceleration I've been trying to tackle this question for a while now, but I'm afraid it's not going anywhere without some outside help. 
So let's say we have some body falling towards a planet without any atmosphere, and let's assume that the initial velocity is not equal to zero. Given these conditions I'm trying to derive a formula which would allow me to determine the exact time of impact or, in fact, time to any given point along the trajectory of the falling body . 
Let's define $r_0$ as initial separation; $r$ as a destination point on the trajectory of the falling body(surface, for example); $v_0$ as initial velocity; $v$ as instantaneous velocity at point $r$;  $t_0$ as initial time; and $T$ as time at point $r$.  
Now acceleration as a function of distance is given by: 
$$\frac{d^2r}{dt^2} = \frac{GM}{r^2}$$
Using the chain rule we arrive at:
$$\ a(r) = \frac{dv}{dt}= \frac{dv}{dr}\frac{dr}{dt}=v\frac{dv}{dr}$$
Therefore
$$a(r)dr= vdv$$
Equivalently:
$$\frac{GM}{r^2}dr=vdv$$
Now in order to get $v$:
$$v = \frac{dr}{dt} = \sqrt{\frac{2G(m_1 + m_2)}{r} - \frac{2G(m_1 + m_2)}{r_0} + v_0^2}$$
$$v = \frac{dr}{dt} =  \sqrt{ \frac{2 G r_0 (m_1 + m_2) - 2 G r (m_1 + m_2)+ v_0^2r_0 r} {r\ r_0}}$$
$dt$ is calculated as: 
$$dt=\sqrt\frac{{r_0 \ r} \ dr} {{ 2 G r_0 (m_1 + m_2) - 2 G r (m_1 + m_2)+(v_0^2rr_0)}}$$
Now I have absolutely no idea how to calculate $dr$, but I hear it has something to do with gravitational potential energy. Some guidance on this matter would be much appreciated. I also have some doubts regarding the next step of my strategy. Let's say I manage to find $dr$, and what then? Intuitively, I feel this is a right path to take but in fact I have no idea how to convert $dr$ and $dt$ into $T$. 

Let me try to help you, but I cannot assure you that my approach is correct. Let $m$ be the mass of falling body and $M$ be the mass of planet. First, we start from the relation of the Newton's second law and Newton's law of universal gravitation.
$$
\begin{align}
ma&=-\frac{GMm}{r^2}\\
\frac{dv}{dt}&=-\frac{GM}{r^2}\\
\frac{dv}{dr}\cdot\frac{dr}{dt}&=-\frac{GM}{r^2}\\
v\ dv&=-\frac{GM}{r^2}\ dr\\
\int_{v_o}^vv\ dv&=-\int_{r_o}^r\frac{GM}{r^2}\ dr\\
\frac{v^2-v_o^2}{2}&=\frac{GM}{r}-\frac{GM}{r_o}\\
v^2&=\frac{2GM}{r}-\frac{2GM}{r_o}+v_o^2\\
v&=\sqrt{\frac{2GM}{r}+v_o^2-\frac{2GM}{r_o}}\\
\frac{dr}{dt}&=\sqrt{\frac{2GM}{r}+v_o^2-\frac{2GM}{r_o}}\\
dt&=\frac{dr}{\sqrt{\frac{2GM}{r}+v_o^2-\frac{2GM}{r_o}}}.
\end{align}
$$
Now, let $a=v_o^2-\dfrac{2GM}{r_o}$ and $b=2GM$, then
$$
\begin{align}
dt&=\frac{dr}{\sqrt{\frac{b}{r}+a}}\\
dt&=\sqrt{\frac{r}{ar+b}}\ dr\\
\int_0^T\ dt&=\int_{r_o}^r\sqrt{\frac{r}{ar+b}}\ dr\\
T&=\int_{r_o}^r\sqrt{\frac{r}{ar+b}}\ dr\\
\end{align}
$$
Let
$$
x^2=\frac{r}{ar+b}\;\quad\Rightarrow\;\quad r= \frac{bx^2}{1-ax^2}\;\quad\Rightarrow\;\quad dr=\frac{2bx}{(1-ax^2)^2}\ dx,
$$
then
$$
T=\int_{r=r_o}^r\frac{2bx^2}{(1-ax^2)^2}\ dx.
$$
The last form integral can be evaluated either using partial fraction expansion, trigonometry substitution, or simply looking the following list of integrals.
$$\\$$

$$\Large\color{blue}{\text{# }\mathbb{Q.E.D.}\text{ #}}$$
