If $R_x$ rotates around the $x$-axis, and $R_y$ rotates around the $y$-axis, and you want to rotate first around $x$, and then around $y$, simply apply $R_y R_x$ to your vector, let's call it $v$.
This is because $R_x v$ rotates $v$ around the $x$-axis, then $R_y(R_x v)$ rotates $R_x v$ around the $y$-axis.
This is called the decomposition into Euler angles. The problem is that the result depends on which order you're planning to do the x-, y-, and z-rotations in.
Here's the drill, though:
Step 1: Convert your rotation to a $3 \times 3$ orthogonal matrix using Rodrigues' formula. This is well-detailed on Wikipedia and elsewhere, so I won't go into it here.
Step 2: Your matrix $M$ has three columns, $u_1, u_2, u_3$. Let's say you want to write
$$
M = R_x R_y R_z
$$
where the $R$s are rotations (by varying amounts) around the three axes. We could then write
$$
R_y^{-1}R_x^{-1}M = R_z
$$
which, since a rotation about $z$ has interesting entries only in its upper left $2 \times 2$ block, means that we need
$$
R_y^{-1}R_x^{-1}u_3 = e_3
$$
where $e_3$ is the vector $\begin{bmatrix}0\\0\\1\end{bmatrix}$.
Let's make that happen. First, we'll find an $x$-rotation that "kills off" the first entry of $u_3$.
Case 1: Both the first and last entries of $u_3$ are 0. In this case, perform no rotation about the $y$ axis, and if the second entry is 1, pick $R_x$ to rotate about $x$ by $90$ degrees; if the second entry is $-1$, pick $R_x$ to rotate about $x$ by $-90$ degrees. With these choices, $R_x^{-1}u_3$ will be $e_3$.
Case 2: At least one of these entries is nonzero. Let's call them $s$ and $c$. Let $\alpha = \text{atan2}(c, s)$, and let $R_y$ be rotation about $y$ by angle $\alpha$. Then $R_y^{-1} u_3$ will have its first entry zero, and its third entry nonzero. Call the second and third entries of $R_y^{-1} u_3$ by the names $a$ and $b$ respectively. Let $\beta = \text{atan2}(b, a)$, and let $R_x$ be rotation about $x$ by angle $\beta$. Then
$$
R_x^{-1}R_y^{-1} u_3 = e_3.
$$
That means that
$$
R_x^{-1}R_y^{-1} M = \begin{bmatrix}
c' & -s' & 0\\
s' & c' & 0\\
0 & 0 & 1
\end{bmatrix},
$$
and we can write $\gamma = \text{atan2}(c', s')$ and let $R_z$ be rotation about $z$ by angle $\gamma$, and we're done.
It's quite possible that I've made a sign error here somewhere: that one of the angles should be the negative of what I wrote. But if you write a program to implement this, it should be pretty obvious where any error is.
Best Answer
A rotation about an arbitrary line $\mathbf u+t\mathbf v$ in $\mathbb R^3$ can be decomposed into a translation to the line, a rotation about an axis through the origin, and a translation back. In homogeneous coordinates, this can be represented by the $4\times4$ matrix $$M=\left[\begin{array}{c|c}I&\mathbf u\\\hline 0&1\end{array}\right]\left[\begin{array}{c|c}R_{\mathbf v,\theta}&0\\\hline 0&1\end{array}\right]\left[\begin{array}{c|c}I&-\mathbf u\\\hline 0&1\end{array}\right]=\left[\begin{array}{c|c}R_{\mathbf v,\theta}&(I-R_{\mathbf v,\theta})\mathbf u\\\hline 0&1\end{array}\right].$$ It’s clear from inspection that this is equivalent to a rotation about a parallel axis that goes through the origin followed by a translation. You have a matrix of the form $$\left[\begin{array}{c|c}R_{\mathbf v,\theta}&\mathbf w\\\hline 0&1\end{array}\right]\tag{*}$$ so recovering the axis of rotation is simply a matter of solving $(I-R_{\mathbf v,\theta})\mathbf u=\mathbf w$ for $\mathbf u$. Unless $\theta$ is an integer multiple of $2\pi$, the matrix $I-R_{\mathbf v,\theta}$ has rank 2, so the solution set—if it exists—will be a line. Note that this equation’s having a solution is a necessary condition for a matrix of the form (*) to represent a rotation. Also, $(I-R_{\mathbf v,\theta})\mathbf u$ is the difference between $\mathbf u$ and its image under $R_{\mathbf v,\theta}$ and so is orthogonal to $\mathbf v$. This gives us another necessary condition: $\mathbf v\cdot\mathbf w=0$.
You can also work directly with the $4\times4$ matrix $M$: $1$ is an eigenvalue and $M$ is the identity on its eigenspace. For most angles $\ker(I-M)$ will be two-dimensional, and so correspond to a line in $\mathbb R^3$. If you use the usual method of row-reduction to find a basis for the kernel, one of the resulting vectors will have $0$ for its last coordinate—it’s a point at infinity. You can take the first three coordinates of this vector as the direction vector $\mathbf v$, while the other vector in the kernel basis can be taken as the fixed vector $\mathbf u$.
The rotation angle can be recovered from either the full matrix $M$ via $\operatorname{tr}M=2+2\cos\theta$, or from the submatrix $R$ via $\operatorname{tr}R=1+2\cos\theta$. These identities come from the facts that the trace of a matrix is equal to the sum of its eigenvalues (with the appropriate number of repetitions), and that the eigenvalues of these matrices are $1$, $e^{i\theta}$ and $e^{-i\theta}$.
Addition: Your idea of using two points and their images to define a pair of planes that intersect along the axis of rotation is sound. It seems like a lot more work than recovering it directly from the matrix to me, though.
As an example, let $\mathbf u=(1,0,-3)^T$, $\mathbf v=(1,-1,-1)^T$ and $\theta = \pi/6$. The rotation submatrix is easily obtained via Rodrigues’ formula and other methods. After multiplying everything out, the resulting matrix is $$M=\begin{bmatrix} {1+\sqrt3\over3}&{-1+\sqrt3\over3}&-\frac13&{-1-\sqrt3\over3} \\ -\frac13&{1+\sqrt3\over3}&{1-\sqrt3\over3}&{4-3\sqrt3\over3} \\ {-1+\sqrt3\over3}&\frac13&{1+\sqrt3\over3}&{-5+2\sqrt3\over3} \\ 0&0&0&1\end{bmatrix}.$$ $\operatorname{tr}M=2+\sqrt3$, so $\cos\theta=\frac{\sqrt3}2$ and so $\theta=\frac\pi6$. Row-reducing $I-M$ produces $$\begin{bmatrix}1&0&1&2\\0&1&-1&-3\\0&0&0&0\\0&0&0&0\end{bmatrix}$$ from which we can read that the rotation axis is $t(-1,1,1)^T+(-2,3,0)^T=(-(t+2),t+3,t)^T$. Setting $t=-3$ gives us the original vector $\mathbf u$. Similarly, row-reducing $\left[\begin{array}{c|c}I-R&\mathbf w\end{array}\right]$ produces $$\left[\begin{array}{ccc|c}1&0&1&-2\\0&1&-1&3\\0&0&0&0\end{array}\right]$$ from which we can see that $\mathbf u=(-(z+2),z+3,z)^T$ as before.