Given the following real input parameters:
$$
\left(x_i,\,y_i\right),
\quad \quad \quad
\left(x_f,\,y_f\right),
\quad \quad \quad
\left(a_x,\,a_y\right),
\quad \quad \quad
v_i > 0
$$
and define the differences:
$$
\Delta x := x_f - x_i\,,
\quad \quad \quad
\Delta y := y_f - y_i
$$
we need to solve the system of equations:
$$
\begin{cases}
\Delta x = v_i\cos\theta\,\Delta t + a_x\,\Delta t^2/2 \\
\Delta y = v_i\sin\theta\,\Delta t + a_y\,\Delta t^2/2 \\
\end{cases}
$$
where solutions $\Delta t > 0$ and $0 \le \theta < 2\pi$ are of interest.
So, rewriting the system of equations:
$$
\begin{cases}
\frac{2\,\Delta x - a_x\,\Delta t^2}{2\,v_i\,\Delta t} = \cos\theta \\
\frac{2\,\Delta y - a_y\,\Delta t^2}{2\,v_i\,\Delta t} = \sin\theta \\
\end{cases}
$$
it follows that:
$$
\left(\frac{2\,\Delta x - a_x\,\Delta t^2}{2\,v_i\,\Delta t}\right)^2 + \left(\frac{2\,\Delta y - a_y\,\Delta t^2}{2\,v_i\,\Delta t}\right)^2 = 1
$$
that is, by reworking:
$$
\left(a_x^2 + a_y^2\right)\Delta t^4 - 4\left(v_i^2 + a_x\,\Delta x + a_y\,\Delta y\right)\Delta t^2 + 4\left(\Delta x^2 + \Delta y^2\right) = 0
$$
where, defined the parameters:
$$
a := a_x^2 + a_y^2\,,
\quad \quad \quad
b := - 4\left(v_i^2 + a_x\,\Delta x + a_y\,\Delta y\right),
\quad \quad \quad
c := 4\left(\Delta x^2 + \Delta y^2\right)
$$
it's clear that there are solutions $\Delta t > 0$ if and only if:
$$
a + c \ne 0
\quad \quad \quad
\text{and}
\quad \quad \quad
b < 0
\quad \quad \quad
\text{and}
\quad \quad \quad
b^2 - 4\,a\,c \ge 0
$$
and in this case it's only one if:
$a = 0$ from which $\Delta t = \sqrt{\frac{c}{-b}}$;
$c = 0$ from which $\Delta t = \sqrt{\frac{-b}{a}}$;
$b^2 - 4\,a\,c = 0$ from which $\Delta t = \sqrt{\frac{-b}{2\,a}}$;
otherwise are two:
$$
\Delta t_1 = \sqrt{\frac{- b - \sqrt{b^2 - 4\,a\,c}}{2\,a}} < \Delta t_2 = \sqrt{\frac{- b + \sqrt{b^2 - 4\,a\,c}}{2\,a}}\,.
$$
In particular, if there are solutions $\Delta t > 0$, the initial velocity vector is trivially determined:
$$
v_{x,i} = v_i\,\cos\theta = \frac{\Delta x}{\Delta t} - \frac{a_x}{2}\,\Delta t\,,
\quad \quad \quad
v_{y,i} = v_i\,\sin\theta = \frac{\Delta y}{\Delta t} - \frac{a_y}{2}\,\Delta t
$$
and consequently also the angle measured with respect to the positive semi-axis of the abscissa:
$$
\theta =
\begin{cases}
\pi + \arctan\left(v_{y,i}/v_{x,i}\right) & \text{if} \; v_{x,i} < 0 \\
3\pi/2 & \text{if} \; v_{x,i} = 0 \, \land \, v_{y,i} < 0 \\
\pi/2 & \text{if} \; v_{x,i} = 0 \, \land \, v_{y,i} > 0 \\
2\pi + \arctan\left(v_{y,i}/v_{x,i}\right) & \text{if} \; v_{x,i} > 0 \, \land v_{y,i} < 0 \\
\arctan\left(v_{y,i}/v_{x,i}\right) & \text{if} \; v_{x,i} > 0 \, \land v_{y,i} \ge 0 \\
\end{cases}
$$
which, if desired, is equivalent to writing:
$$
\theta =
\begin{cases}
2\pi + \text{arctan2}\left(v_{y,i},\,v_{x,\,i}\right) & \text{if} \; v_{y,i} < 0 \\
\text{arctan2}\left(v_{y,i},\,v_{x,\,i}\right) & \text{if} \; v_{y,i} \ge 0 \\
\end{cases}
$$
where $\text{arctan2}(y,\,x)$ associates an angle in the interval $(-\pi,\,\pi]$ for each point $(x,\,y) \ne (0,\,0)$.
Since this is a uniformly accelerated motion, both in the horizontal and vertical directions, the parabolic trajectories will generally be rototranslated in the plane. For example, assuming:
$$
\left(x_i,\,y_i\right) = (1,\,1),
\quad \quad \quad
\left(x_f,\,y_f\right) = (3,\,2),
\quad \quad \quad
\left(a_x,\,a_y\right) = (1,\,-1),
\quad \quad \quad
v_i = 2
$$
the two trajectories solution of the problem are:
where blue is for $\Delta t_1 \approx 1.06$ and $\theta_1 \approx 47.4°$, while red is for $\Delta t_2 \approx 2.98$ and $\theta_2 \approx 114°$.
Finally, it isn't difficult to prove that the area enclosed by the two parabolic arcs is equal to:
$$
\mathcal{A} = \frac{\left|a_x\,\Delta y - a_y\,\Delta x\right|\left(\Delta t_2^2 - \Delta t_1^2\right)}{12} = \frac{\sqrt{15}}{2} \approx 1.94\,.
$$
I hope it's clear enough, good luck! ^_^
