Using Kjetil's answer answer, with process91's comment, we arrive at the following procedure.

## Derivation

We are given two unit column vectors, $A$ and $B$ ($\|A\|=1$ and $\|B\|=1$). The $\|\circ\|$ denotes the L-2 norm of $\circ$.

First, note that the rotation from $A$ to $B$ is just a 2D rotation on a plane with the normal $A \times B$. A 2D rotation by an angle $\theta$ is given by the following augmented matrix:
$$G=\begin{pmatrix}
\cos\theta & -\sin\theta & 0 \\
\sin\theta & \cos\theta & 0 \\
0 & 0 & 1
\end{pmatrix}.$$

Of course we don't want to actually *compute* any trig functions. Given our unit vectors, we note that $\cos\theta=A\cdot B$, and $\sin\theta=||A\times B||$. Thus
$$G=\begin{pmatrix}
A\cdot B & -\|A\times B\| & 0 \\
\|A\times B\| & A\cdot B & 0 \\
0 & 0 & 1\end{pmatrix}.$$

This matrix represents the rotation from $A$ to $B$ in the base consisting of the following column vectors:

normalized vector projection of $B$ onto $A$: $$u={(A\cdot B)A \over \|(A\cdot B)A\|}=A$$

normalized vector rejection of $B$ onto $A$: $$v={B-(A\cdot B)A \over \|B- (A\cdot B)A\|}$$

the cross product of $B$ and $A$: $$w=B \times A$$

Those vectors are all orthogonal, and form an orthogonal basis. This is the detail that Kjetil had missed in his answer. You could also normalize $w$ and get an orthonormal basis, if you needed one, but it doesn't seem necessary.

The basis change matrix for this basis is:
$$F=\begin{pmatrix}u & v & w \end{pmatrix}^{-1}=\begin{pmatrix} A & {B-(A\cdot B)A \over \|B- (A\cdot B)A\|} & B \times A\end{pmatrix}^{-1}$$

Thus, in the original base, the rotation from $A$ to $B$ can be expressed as right-multiplication of a vector by the following matrix: $$U=F^{-1}G F.$$

One can easily show that $U A = B$, and that $\|U\|_2=1$. Also, $U$ is the same as the $R$ matrix from Rik's answer.

## 2D Case

For the 2D case, given $A=\left(x_1,y_1,0\right)$ and $B=\left(x_2,y_2,0\right)$, the matrix $G$ is the forward transformation matrix itself, and we can simplify it further. We note
$$\begin{aligned}
\cos\theta &= A\cdot B = x_1x_2+y_1y_2 \\
\sin\theta &= \| A\times B\| = x_1y_2-x_2y_1
\end{aligned}$$

Finally,
$$U\equiv G=\begin{pmatrix}
x_1x_2+y_1y_2 & -(x_1y_2-x_2y_1) \\
x_1y_2-x_2y_1 & x_1x_2+y_1y_2
\end{pmatrix}$$
and
$$U^{-1}\equiv G^{-1}=\begin{pmatrix}
x_1x_2+y_1y_2 & x_1y_2-x_2y_1 \\
-(x_1y_2-x_2y_1) & x_1x_2+y_1y_2
\end{pmatrix}$$

## Octave/Matlab Implementation

The basic implementation is very simple. You could improve it by factoring out the common expressions of `dot(A,B)`

and `cross(B,A)`

. Also note that $||A\times B||=||B\times A||$.

```
GG = @(A,B) [ dot(A,B) -norm(cross(A,B)) 0;\
norm(cross(A,B)) dot(A,B) 0;\
0 0 1];
FFi = @(A,B) [ A (B-dot(A,B)*A)/norm(B-dot(A,B)*A) cross(B,A) ];
UU = @(Fi,G) Fi*G*inv(Fi);
```

Testing:

```
> a=[1 0 0]'; b=[0 1 0]';
> U = UU(FFi(a,b), GG(a,b));
> norm(U) % is it length-preserving?
ans = 1
> norm(b-U*a) % does it rotate a onto b?
ans = 0
> U
U =
0 -1 0
1 0 0
0 0 1
```

Now with random vectors:

```
> vu = @(v) v/norm(v);
> ru = @() vu(rand(3,1));
> a = ru()
a =
0.043477
0.036412
0.998391
> b = ru()
b =
0.60958
0.73540
0.29597
> U = UU(FFi(a,b), GG(a,b));
> norm(U)
ans = 1
> norm(b-U*a)
ans = 2.2888e-16
> U
U =
0.73680 -0.32931 0.59049
-0.30976 0.61190 0.72776
-0.60098 -0.71912 0.34884
```

## Implementation of Rik's Answer

It is computationally a bit more efficient to use Rik's answer. This is also an Octave/MatLab implementation.

```
ssc = @(v) [0 -v(3) v(2); v(3) 0 -v(1); -v(2) v(1) 0]
RU = @(A,B) eye(3) + ssc(cross(A,B)) + \
ssc(cross(A,B))^2*(1-dot(A,B))/(norm(cross(A,B))^2)
```

The results produced are same as above, with slightly smaller numerical errors since there are less operations being done.