For fun, I am trying to find the matrix which corresponds to a rotation of angle $\theta$ around an unit axis vector $\hat{n}$ using calculus. Here are the intuitive properties we begin with:
$$ \begin{align} &1. M_{\theta} M_{\phi} = M_{\phi} M_{\theta} \\ &2. M_{\phi} n = n \\ &3. M_{\phi + \theta} = M_{\phi} M_{\theta}\end{align}$$
I begin with computing the derivative of rotation matrix with $\phi$:
$$ \frac{d}{d \phi} M_{\phi} = \lim_{d \phi \to 0} \frac{M_{\phi + d \phi} - M_{\phi} } {d \phi} = \lim_{d \phi \to 0} \frac{M_{\phi} ( M_{d \phi} - I) }{d \phi}= M_{\phi} \lim_{d \phi \to 0} (\frac{M_{d \phi} - I}{d \phi}) $$
I need to figure out the last derivative... a quick thought was that it reminded me of the steps involved in taking derivative of exponential function where in the intermediate step we have to deal with the limit $\lim_{d \phi \to 0} \frac{e^{ d \phi} - 1}{d \phi}$.
However the observation doesn't help much to do computations. I thought of that we can think instead of how matrix acts on vectors to compute the derivative. Consider a vector $\vec{v}$ it can be decomposed into component perpendicular and parallel to $n$. $$ v= v_n \hat{n} + (v-v_n \hat{n} ) $$
Now, $$(M_{d \phi} - I) v= (M_{d \phi} - I)( v - v_{n} \hat{n})$$
By geometry, we can find the above expression must be equal to $ n \times (v-v_{n} \hat{n} ) = n \times v $. Picture:
I now have:
$$ \frac{d}{d \phi} (M_{\phi} v)= M_{\phi} \lim_{d \phi \to 0} (\frac{\left[ M_{d \phi} - I \right] v}{d \phi}) = M_{\phi} (n \times v)$$
Now I need to cancel the $v$ on both sides. One way would be to associate $n$ to it's cross product matrix, let that be $B$:
$$ \frac{d}{d \phi} (M_{\phi} v) =( M_{\phi}) B v$$
$$ \frac{d}{d \phi} M_{\phi} = M_{\phi} B$$
A possible continuation:
The solution to $ \frac{d}{d \phi} M_{\phi} = M_{\phi} B$, we can guess to be:
$$ M_{\phi} = e^{\phi B}$$
Since, the solution to $ \frac{dy}{dx} = ya $ is $y=e^{xa}$. By the power series, we have:
$$ M_{\phi} = I + \phi B + \frac{(\phi B)^2}{2!}....$$
The idea behind powers of $B$ is that we cross product as many times as the power. Eg: $B^3v = n \times ( n \times ( n \times v ))$.
Note: the above is actually computable through cross products without knowing how to convert cross product into matrix format. We can now conver the formula to look more like eulers formula. Consider an orthonormal basis $p,q,n$ satisfying:
$$ p \times q = n $$
$$ n \times p = q$$
$$ q \times n = p $$
Now, the consider a general vector $v$ which $M_{\phi}$ can act on, it can be written as: $ v = v_n n + v_p p + v_q q = v_n n + v'$, we now observe the pattern as powers of $B$ acts on it:
$$B^0 v= v = v_n n + v'$$
$$B^1 v = v_p q - v_q p=B^1 v' $$
$$ B^2 v = -v_p p- v_q q = -v'$$
$M_{\phi}$ action on a vectro: $$M_{\phi}v = (I + \phi B + \frac{\phi^2 B^2}{2!} + \frac{\phi^3 B^3}{3!}...)v = (I + \phi B + \frac{\phi^2 B^2}{2!} + \frac{\phi^3 B^3}{3!}...)(v_n n + v') \\= v_n n+ (I + \phi B + \frac{\phi^2 B^2}{2!} + \frac{\phi^3 B^3}{3!}...)( v') \\=v_n n + ( I+ \phi B - \frac{\phi^2}{2!} I - \frac{\phi^3 B}{3!}...) v' \\ = v_n n + (\cos \phi +B \sin \phi)v'\\=v+B^2 v - ( \cos \phi + B \sin \phi)B^2 v$$
Cancelling $v$:
$$M_{\phi}= I+B^2 -(\cos \phi + B \sin \phi) B^2= I +( 1- \cos \phi - B \sin \phi)B^2$$
Questions:
- Is my final answer correct?
- Has this been done before?
- Is it possible to make my arguments rigorous?

I think the direct approach using standard vector algebra is much better in this case than calculus. You can derive the rotation matrix very easily as follows:
Let the axis of rotation be the unit vector $a$, and let the vector you want to rotate be $v$.
Decompose $v$ into $v_1$ along the axis and $v_2$ perpendicular to the axis, then
$v_1 = (v \cdot a) a = a a^T v $ and $v_2 = v - v_1 = (I - a a^T) v $
Construct the vector that has the same magnitude as $v_2$ but orthogonal to both the axis and $v_2$. This vector is $v_3 = a \times v_2 = a \times (I - a a^T) v = a \times v$
Now the rotated vector is a combinaton of $v_1, v_2$ and $v_3$ as follows:
$v' = R v = v_1 + \cos \theta v_2 + \sin \theta v_3 $
Substitute what $v_1, v_2, v_3$ are and you're almost done. Hence,
$v' = R v = a a^T v + \cos \theta (I - a a ^T ) v + \sin \theta (a \times v) \hspace{20pt} (1) $
This last equation is known as the Rodrigues' rotation formula for a vector $v$ about an axis $a$. To extract the rotation matrix we need to express $a \times v$ in matrix form. Writing the components of $a \times v$ in terms of the components of $a$ and $v$, one deduces that
$ a \times v = S_a v \hspace{20pt} (2) $
where
$S_a = \begin{bmatrix} 0 && - a_z && a_y \\ a_z && 0 && - a_x \\ -a_y && a_x && 0 \end{bmatrix} $
Using this last equation $(2)$ into $(1)$ and factoring $v$ out, we get the rotation matrix as
$ R = a a ^T + ( I - a a ^T ) \cos \theta + S_a \sin \theta \hspace{20pt}(4) $
That's all to it.
Now, let's check if the matrix derived by the OP matches what I got here.
The final matrix obtained by the OP is
$ M_\phi = I + ( (1 - \cos \phi) I - B \sin \phi ) B^2 $
Now, $B^2 v = n \times (n \times v) = n (n^T v) - v (n^T n) = n n^T v - v = (n n^T - I ) v$
Replacing $B^2 $ with $(n n^T - I)$ in the above formula gives
$ M_\phi = I + ( (1 - \cos \phi) I - B \sin \phi ) (n n^T - I ) $
This expands into
$ M_\phi = I + (1 - \cos \phi) (n n^T - I) + B \sin \phi $
(Because $B n = 0 $ )
Simplifing further,
$ M_\phi = I + n n^T - I + \cos \phi (I - n n^T ) + B \sin \phi $
And this equal to
$ M_\phi = n n^T + \cos \phi (I - n n^T ) + B \sin \phi \hspace{20pt}(*)$
Comparing this with the result I derived earlier ( equation $(4)$ ) , we find that $(*)$ is identical to $(4)$ once we make the substitutions $n = a$ , $B = S_a$ and $\phi = \theta$.