There are a number of ways to do so. Using quaternions has already been mentioned so I will go with (4x4) transformation matrices.
Let a transformation matrix $T$ be
$$
T = \begin{bmatrix}R & p\\ 0 & 1\end{bmatrix},
$$
where $R$ is a rotation matrix and $p$ is a position vector. I will write $T = (R, p)$ for short.
You want to generate a path connecting $T_1= (R_1, p_1)$ and $T_2= (R_2, p_2)$. The problem then boils down to generating a rotational path and a translational path, separately. For translation, a path can be generated pretty easily using any method like polynomial interpolation or using splines. For rotation, it's a bit more complicated but in the end it can also be done using methods like polynomial interpolation. So I will talk a bit about rotation here.
The thing is that any rotation matrix can be written as
$$R = e^{[r]},$$
where $r \in \mathbf{R}^3$, $[\cdot]$ maps a vector to a unique skew symmetric matrix, and $e^{[r]}$ is the matrix exponential
$$
e^{[r]} = \sum_{k = 0}^{\infty}\frac{1}{k!}[r]^k.
$$
(Of course, there is a formula for that without evaluating that infinite sum.) So actually a rotational path may be generated as
$$R(t) = R_0e^{[r(t)]},$$
where $r(t)$ is a polynomial vector and $R_0$ is some constant rotation matrix. Then you just need to find $r(t)$ such that $R_0e^{[r(0)]} = R_1$ and $R_0e^{[r(T)]} = R_2$, so as to meet your boundary conditions. For more details, you may have a look at this paper, for example.
