There is a direct formula for the rotation matrix for an arbitrary axis and angle. Given a unit vector $a = (a_x, a_y, a_z)$ and angle $\theta$, the matrix can be constructed as follows (derivation from Wikipedia):
First build a matrix $C$ from the components of $a$ according to the following formula:
$$ C =
\begin{bmatrix}
0 & -a_z & a_y \\
a_z & 0 & -a_x \\
-a_y & a_x & 0
\end{bmatrix} $$
Then you can construct the rotation matrix from $C$ as follows:
$$ R_a(\theta) = I + C \sin \theta + C^2 (1 - \cos \theta) $$
where $I$ is the 3×3 identity matrix. This assumes you're using a column-vector convention; if you're using row vectors instead, transpose $C$.
By the way, this formula has a neat geometric interpretation. The matrix $C$ has the effect of calculating cross products with $a$. In other words, for any vector $v$, $Cv = a \times v$. If you imagine rotating some vector $v$ about $a$, the tip of $v$ traces out a circle perpendicular to $a$; $Cv$ is a vector tangent to that circle, and $C^2 v$ (which equals $a \times (a \times v)$) is a vector normal to the circle. So with those two vectors, you have a basis for the plane of the circle; then you can apply a variation of the usual sine/cosine formula for points on a circle.
