Utilizing maximum acceleration $a$ for displacement $d$ with initial velocity $v_0$ and final velocity $v_1$ Problem
My goal is to move an object from point a to b (displacement $d$) as fast as possible utilizing the maximum available acceleration $a_{max}$, taking into account the initial velocity $v_0$ and final velocity $v_1$.
Elaboration
I'm working on spaceship physics which means there is no friction or drag. The direction of the object is not necessarily parallel to the direction of it's velocity, neither are the initial and final velocities necessarily parallel.

Solving for one dimension
I started off by solving this in one dimension using the formulas for constant acceleration, specifically:
$$
d = \frac{v_1^2-v_0^2}{2a}
$$
As my intuition tells me there should be a displacement $d_0$ where constant acceleration $a_{max}$ is performed (where the object reaches its maximum velocity $v_{max}$) followed by a displacement $d_1$ where constant deceleration $-a_{max}$ is performed. From the above equation, and from:
$$
d = d_0 + d_1
$$
We find that:
$$
d = \frac{v_{max}^2-v_0^2}{2a_{max}} + \frac{v_1^2-v_{max}^2}{-2a_{max}}
$$
Solving for $v_{max}$ gives:
$$
v_{max}^2 = \frac{2a_{max}d+v_1^2+v_0^2}2
$$
Using this information I can accelerate the object until $v_{max}$ is reached and then decelerate it.

Solving for 3 dimensions
Doing this in 3-dimensions is a bit more difficult. I'm assuming that $||\vec{a}||$ should be constant and equal the maximum available acceleration $a_{max}$.
My first attempt here was to calculate $v_{max}$ (like above) for each axis seperately. The problem though is to find a proper acceleration distribution along each axis when the displacement changes (having a moving target), so that the object for example does not go into orbit.
In another attempt I did the following spring/damper like calculation to find $\vec{a}$:
$$
\vec{s} = \vec{v_1}-\vec{v_0}+\vec{d}
$$
$$
\vec{a} = \frac{\vec{s}}{||\vec{s}||} * a_{max}
$$
This works fine when accelerating but deceleration starts off too early.
I have been trying to search the web for this information but to no avail, any suggestions?
 A: I ended up using a solution by Christer Swahn from his blog: http://mmoarch.blogspot.se/2012/05/computing-space-travel.html
The trajectory is approximated and then optimised using the bisection method. A near perfect result is usually reached after 10 iterations.
See full solution (Java implementation) at: http://mmoarch.blogspot.se/2012/05/computing-space-travel-implementation.html
A: Problem. On object of mass $m$
  is located at $\vec{r}_{A}=\begin{pmatrix}0 & 0 & 0\end{pmatrix}$
  with velocity $\vec{v}_{A}=\begin{pmatrix}v_{A} & 0 & 0\end{pmatrix}$
  and needs to reach point $\vec{r}_{B}=\begin{pmatrix}d_{x} & d_{y} & 0\end{pmatrix}$
  with velocity $\vec{v}_{B}=\begin{pmatrix}v_{B}\cos\theta & v_{B}\sin\theta & 0\end{pmatrix}$
  in the minimum amount of time when acceleration magnitude is limited to $\left|\vec{a}\right|\leq g$
 .
Let us proceed with a general solution and try to optimize it in the end. I propose a five step process (inspired from race car driving). I change the notation to keep track of the steps better. Also the direction vector at A
  is $\vec{e}_{A}=\begin{pmatrix}1 & 0 & 0\end{pmatrix}$
  and at B
  is $\vec{e}_{B}=\begin{pmatrix}\cos\theta & \sin\theta & 0\end{pmatrix}$
 . Also the normal vectors are $\vec{n}_{A}=\begin{pmatrix}0 & 1 & 0\end{pmatrix}$
  and $\vec{n}_{B}=\begin{pmatrix}-\sin\theta & \cos\theta & 0\end{pmatrix}$
 , where $\theta$
 is the required turn angle between A
  and B
 . This can be found by $\cos\theta=\pm\frac{\vec{v}_{A}\cdot\vec{v}_{B}}{\left|\vec{v}_{A}\right|\left|\vec{v}_{B}\right|}$
 . To find the sign use ${\rm sign}\left(\hat{k}\cdot\left(\vec{v}_{A}\times\vec{v}_{B}\right)\right)$
  to see if you need to turn left or right.


*

*Accelerate to a maximum speed $v_{1}$
on a straight line from the starting location $\vec{r}_{0}=\vec{r}_{A}$
and speed $\vec{v}_{0}=v_{A}\vec{e}_{A}$
. The final speed is $\vec{v}_{1}=v_{1}\vec{e}_{A}$
, position $\vec{r}_{1}=\vec{r}_{A}+\frac{v_{1}^{2}-v_{A}^{2}}{2g}\vec{e}_{A}$
, and time $t_{1}=\frac{v_{1}-v_{A}}{g}$
.

*Decelerate to a desired cornering speed $v_{C}$
on a straight line. The final speed is $\vec{v}_{2}=v_{C}\vec{e}_{A}$
, position $\vec{r}_{2}=\vec{r}_{1}+\frac{v_{1}^{2}-v_{C}^{2}}{2g}\vec{e}_{A}$
, and time $t_{2}=t_{1}+\frac{v_{1}-v_{C}}{g}$
. 

*Turn to align your speed with the final velocity $\vec{v}_{B}$
and position along the straight line defined by the final location $\vec{r}_{B}$
and velocity direction. The turn radius is $r=\frac{v_{C}^{2}}{g}{\rm sign}(\theta)$
. The final speed is $\vec{v}_{3}=v_{C}\vec{e}_{B}$
, position is $\vec{r}_{3}=\vec{r}_{2}+\left(\vec{n}_{A}-\vec{n}_{B}\right)\frac{v_{C}^{2}}{g}{\rm sign}(\theta)$
, and time $t_{3}=t_{2}+\frac{v_{C}}{g}\theta$
.

*Accelerate to a maximum speed $v_{4}$
on a straight line. The final speed is $\vec{v}_{4}=v_{4}\vec{e}_{B}$
, position $\vec{r}_{4}=\vec{r}_{3}+\frac{v_{4}^{2}-v_{C}^{2}}{2g}\vec{e}_{B}$
, and time $t_{4}=t_{3}+\frac{v_{4}-v_{C}}{g}$
.

*Decelerate to the final location $\vec{r}_{B}$
and speed $\vec{v}_{B}=v_{B}\vec{e}_{B}$
on a staight line. The final speed is $\vec{v}_{B}$
, position $\vec{r}_{B}=\vec{r}_{4}+\frac{v_{4}^{2}-v_{B}^{2}}{2g}\vec{e}_{B}$
, and time $t_{B}=t_{4}+\frac{v_{4}-v_{B}}{g}$
.

The total time is $$t=\frac{1}{g}\left(2v_{1}+2v_{4}+\left(\theta-2\right)v_{C}\right)-\frac{1}{g}\left(v_{A}+v_{B}\right)$$
 where $v_{1}$ , $v_{4}$  and $v_{C}$
  are unknown. The final position is $$d_{x} = \frac{2v_{1}^{2}-v_{A}^{2}-v_{C}^{2}}{2g}+\frac{2v_{4}^{2}-v_{B}^{2}-v_{C}^{2}}{2g}\cos\theta+\frac{v_{C}^{2}}{g}\left|\sin\theta\right|$$
 $$d_{y}=\frac{2v_{4}^{2}-v_{B}^{2}-v_{C}^{2}}{2g}\sin\theta+\frac{v_{C}^{2}}{g}\left(1-\cos\theta\right){\rm sign}\theta$$
The final position is used to find speeds $v_{1}$
  and $v_{4}$
 $$v_{1}=\sqrt{\frac{v_{A}^{2}}{2}+\frac{v_{C}^{2}}{2}+\frac{g\left(d_{x}\sin\theta-d_{y}\cos\theta\right)}{\sin\theta}-v_{C}^{2}\left|\tan\frac{\theta}{2}\right|}$$
 $$v_{4}=\sqrt{\frac{v_{B}^{2}}{2}+\frac{v_{C}^{2}}{2}+\frac{g\; d_{y}}{\sin\theta}-v_{C}^{2}\left|\tan\frac{\theta}{2}\right|}$$
In the end the time is a function of $v_{C}$
  only but it is not easy to optimize this function. Consider solving the following
$$ \frac{{\rm d} t}{{\rm d} v_C} = \frac{1}{g} \left( 2 \frac{{\rm d} v_1}{{\rm d} v_C} + 2 \frac{{\rm d} v_4}{{\rm d} v_C} + \theta-2 \right) = 0 $$
