Given 2 constant acceleration model with initial position in 2D, how can I compute the when will the two object be closest to each other? The motion model is defined by
Initial position $\vec{x}$ in 2d space
Initial timestamp $t$
Initial velocity $\vec{v}$ of $x,y$ direction
a constant acceleration $\vec{a}$ of $x,y$ direction
How can I find the timestamp $t$ that, two object $(\vec{x}_0, \vec{v}_0, \vec{a}_0, t_0)$ and $(\vec{x}_1, \vec{v}_1, \vec{a}_1, t_1)$ are/were closest to each other? 
What I have try:
Let $p_0$ and $p_1$ be the point that they are/were closest to each other at some timestamp $t$ which can be calculated from $p_0=x_0+v_0(t_0-t)+\frac{1}{2}a_0(t_0-t)^2$
I tried to solve $t$ by throwing this into matlab/python sovler $\frac{d|p_0-p_1|_{2}}{dt}=0$. but it give me a huge mess and I am wondering if there is a simpler way to do it?
Edit Note: the two point might not collide
 A: I think the formula you want for predicting positions is
$$\vec p_0 = \vec x_0 + \vec v_0(t - t_0) + \tfrac12 \vec a_0(t - t_0)^2.$$
(The object tends to move in the direction of $\vec v_0$ as time increases.)
But you can simplify things by rewriting the conditions so the initial timestamp is the same.
For example, you can make the initial timestamp $t_0$ in both cases by replacing $\vec x_1$ and $\vec v_1$:
\begin{align}
\vec x'_1 &=  \vec x_1 + \vec v_1(t_0-t_1) + \tfrac12 \vec a_1(t_0-t_1)^2, \\
\vec v'_1 &=  \vec v_1 + \vec a_1(t_0-t_1).
\end{align}
So now you can describe the second body as 
$(\vec x'_1, \vec v'_1, \vec a_1, t_0).$
You can change the problem to one of relative position/velocity of the two objects.
Let 
\begin{align}
\vec x_2 &= \vec x'_1 - \vec x_0,\\
\vec v_2 &= \vec v'_1 - \vec v_0,\\ 
\vec a_2 &= \vec a_1 - \vec a_0.
\end{align}
Then the relative position of the two objects at any time $t$ is
$$\vec p_2(t) = \vec x_2 + \vec v_2(t - t_0) + \frac12 \vec a_2(t - t_0)^2.$$
You want to minimize the length of the vector $\vec p_2(t).$
This amounts to finding the point on a parabola closest to the point $(0,0).$
You have the point on the parabola as a function of time, the slope of the normal also obtainable as a function of time, and you need the line with that slope through that point to pass through $(0,0).$
