Let the point $p$ be given in cartesian coordinates by $(\sin\theta\cos\phi, \sin\theta\sin\phi, \cos\phi)$, so that $\theta$ and $\phi$ are the usual spherical coordinates on $S^2$. Then the matrix that rotates $n=(1,0,0)$ to $p$, and also preserves the plane passing through both of them and the origin, is given by
$M=\left(\begin{array}{ccc}\sin^2\phi + \cos\theta\cos^2\phi & (-1+\cos\theta)\sin\phi\cos\phi & \sin\theta\cos\phi\\ (-1+\cos\theta)\sin\phi\cos\phi & \cos^2\phi+\cos\theta\sin^2\phi & \sin\theta\sin\phi\\ -\sin\theta\cos\phi & -\sin\theta\sin\phi & \cos\theta \end{array}\right).$
It's straightforward but tedious to check that this is indeed an orthogonal matrix with determinant $1$ that sends $n$ to $p$ and fixes $u=n\times p=(-\sin\theta\sin\phi,\sin\theta\cos\phi,0)$, and hence the desired transformation.
I produced this matrix by considering that, for each fixed $\phi$, what you desire is a one-parameter subgroup of $SO(3)$: the group of rotations that fix $u$. It therefore is the family of exponentials of an infinitesimal transformation, i.e. an element of the Lie algebra $\frak{so}\mathrm{(3)}$. For example, a tiny rotation tilting the $z$-axis toward the $x$-axis differs from the identity by $\left(\begin{array}{ccc} & & \epsilon\\ & & \\ -\epsilon & & \end{array}\right)$, and
$\exp\left(\begin{array}{ccc} & & \theta\\ & & \\ -\theta & & \end{array}\right) = \left(\begin{array}{ccc} \cos\theta & & \sin\theta\\ & 1 & \\ -\sin\theta & & \cos\theta \end{array}\right).$
Similarly, the Lie algebra element corresponding to an infinitesimal rotation of the $z$-axis in the direction of the $y$-axis is given by $\left(\begin{array}{ccc} & & \\ & & \epsilon\\ & -\epsilon& \end{array}\right)$, and
$\exp\left(\begin{array}{ccc} & & \\ & & \theta\\ & -\theta & \end{array}\right) = \left(\begin{array}{ccc} 1 & & \\ & \cos\theta & \sin\theta\\ & -\sin\theta & \cos\theta \end{array}\right).$
We want the Lie algebra element which has components $\cos\phi$ toward the $x$-axis and $\sin\phi$ toward the $y$-axis; thus we desire the exponential
$M=\exp\left(\begin{array}{ccc} & & \theta\cos\phi\\ & & \theta\sin\phi\\ -\theta\cos\phi & -\theta\sin\phi & \end{array}\right)=\exp(A).$
This matrix $A$ is diagonalizable, so its exponential $M$ can be computed by hand, which is how I arrived at my answer above.