I have $k$ vectors $v_i\in\mathbb{R}^n$, mutually orthogonal. I would now like to rotate them in the $k$-dimensional subspace spanned by the $v_i$ such that $v_0$ ends up at the given target vector $w\in span(v_i)$.
How do I find the orthonormal rotation matrix $R$?
Short answer: Gram-Schmidt. Longer answer:
Normalize your input vectors, i.e., replace $v_i$ with $v_i / \|v_i\|$ for each $i$.
Create a set of $n+k$ vectors $v_1, v_2, \ldots, v_k, e_1, e_2, \ldots, e_n$. Call these $u_1, \ldots, u_{n+k}$.
Perform the Gram-Schmidt orthogonalization process on each of these. Gram-Schmidt, on the $i$th vector, produces an intermediate vector $b_i$, which is then normalized to produce the output vector $w_i$. $b_i$ is defined recursively as follows:
a. $b_1 = u_1$.
b. For $i > 1$, $b_i = u_i - (u_i \cdot w_1) w_1 - (u_1 \cdot w_2) w_2 - ... - (u_1 \cdot w_{i-1}) w_{i-1}.
If it turns out that $b_i$ is zero after this step, simply discard it and move on to the next $u$.
When you're done, you'll have exactly $n$ non-discarded vectors, and they'll be an orthonormal basis for $\Bbb R^n$, and the first $k$ will span the subspace spanned by the $u_i$s. Put them in a matrix $M$, each as a column of $M$.
Then the matrix $R$ that you're looking for is (almost) $M^t$. The transformation $x \mapsto M^t x$ sends $v_0$ to $e_1$, $v_1$ to $e_2$, and so on, up to $v_k$ being sent to $e_{k+1}$.
Now all we need is an orthogonal transformation that sends $e_1$ to $w$ and the span of $e_1, \ldots, e_{k+1}$ to the span of the $v$ vectors. To do this...
Apply the same modified Gram-Schmidt to the set $\{w, v_1, v_2, \ldots, v_k, e_1, e_2, \ldots, e_n\}$ to produce $n$ vectors $c_1 = w, \ldots, c_n$ that you place as columns of a matrix $C$. Note that $Ce_1 = w$. Also notice that $Ce_i$, for $i = 1, \ldots, k+1$ is in the span of the $v$-vectors.
The matrix $CM^t$ is then the one that you need.
N.B.:
$e_i$ here denotes the $i$th standard basis vector: a column vector of all $0$s except for a $1$ in the $i$th spot.
This approach is not particularly numerically stable; if you discard $b_i$ whenever it's "too short" (and the precise notion of "too short" depends on $n$ in a way I'm not willing to work out), you'll still get a nice orthonormal basis with far better numerical stability.