First off, there is no simple rule (at least as far as I know) to represent a rotation as a cascade of single axis rotations. One mechanism is to compute a rotation matrix and then compute the single axis rotation angles. I should also mention that there is not a unique solution to the problem you are posing, since there are different cascades of rotations that will yield the same final position.
Let's first build two references frames, $A$ and $B$.
$A$ will denote the standard Euclidian basis. $B$ will denote the reference frame which results from applying the desired rotation. I will use subscript notation when referring to a vector in a specific reference frame, so $v_A$ is a vector whose coordinates are defined with respect to the standard Euclidian basis, and $v_B$ is defined w.r.t. $B$.
All we know so far about $B$ is that the vector $(0,0,1)_B = \frac{1}{\sqrt{77}}(8, -3, 2)_A$. We have one degree of freedom to select how $(0,1,0)_B$ is defined in reference frame $A$. We only require that it be orthogonal to $\frac{1}{\sqrt{77}}(8, -3, 2)_A$. So we would like to find a unit vector such that
$$8x-3y+2z=0$$
One solution is $\frac{1}{\sqrt{6}}(1,2,-1)_A$, which you can verify. This forces our hand for how the vector $(1,0,0)_B$ is represented w.r.t. $A$, hence
$$(1,0,0)_B=\frac{1}{\sqrt{462}}(-1,10,19)_A$$
So the rotation matrix from $B$ to $A$ is
$$R_B^A=\left( \begin{array}{ccc}
\frac{-1}{\sqrt{462}} & \frac{1}{\sqrt{6}} & \frac{8}{\sqrt{77}} \\
\frac{10}{\sqrt{462}} & \frac{2}{\sqrt{6}} & \frac{-3}{\sqrt{77}} \\
\frac{19}{\sqrt{462}} & \frac{-1}{\sqrt{6}} & \frac{2}{\sqrt{77}} \end{array} \right)$$
Given any vector with coordinates defined in the Euclidian basis (such as a point in the $xy$-plane), multiplying by the inverse this matrix will yield the representation in our new reference frame. If you want to search for the three axis rotations, I refer you to Wolfram MathWorld's page on Euler Angles.
Hope this helps! Let me know if anything is unclear or if I you think you see any mistakes.
Best Answer
What I understand is that the rotation should rotate the vector $u$ to a vector $u'$ pointing in the direction of the positive $z$ axis.
In the usual interpretation of what a vector is in $\mathbb R^3,$ the vector $u$ does not actually have a particular start or end point; it merely gives a direction and distance; but you can perfectly well specify the vector $u$ by specifying two points, $\newcommand{\startp}{\mathit{start}\,}\startp$ and $\newcommand{\endp}{\mathit{end}\,}\endp$, and taking the differences in those points' coordinates to get the coordinates of $u.$ Moreover, once you have a rotation that makes the image $u'$ point in the positive $z$ direction, the same rotation will take $\startp$ to a point $\startp'$ and take $\endp$ to a point $\endp'$ such that if you name the $x,y,z$ coordinates of these points by using the subscripts $x,y,z$ then \begin{align} \endp_x' &= \startp_x', \\ \endp_y' &= \startp_y', \\ \endp_z' &= \startp_z' + \lVert u\rVert > \startp_z'. \\ \end{align} This seems to match your sketch.
Assuming that rotation matrices are applied on the left side of a vector (the usual convention), the function that you have chosen to use, $$ M_{\mathit{rot}} = M_{Z} M_{Y}M_{X},$$ performs the $Z$ rotation last. This rotation preserves the angles that vectors make with the $z$ axis, so if we are to get $u'$ parallel to the $z$ axis we are not going to do it with this last step; we must do it with the rotations around the $x$ and $y$ axes instead.
The rotation around the $y$ axis will not be able to "correct" any difference in the $y$ coordinates of the start and end points, so we need the first rotation, the one around the $x$ axis, to rotate the points to the same $y$ coordinate. We do this by rotating $u$ to a vector parallel to the $x,z$ plane; that is, eliminate its $y$ coordinate.
You can consider a rotation around the $x$ axis to leave all $x$ coordinates unchanged but to rotate the projections of points on the $y,z$ plane. The projection of $u$ on the $y,z$ plane is $(u_y, u_z)$; we want to rotate this projected vector to the vector $\left(0, \sqrt{u_y^2 + u_z^2}\right),$ which has the same length but whose direction is the positive $z$ direction.
You can find the angle that the vector $(u_y, u_z)$ makes with the $z$ axis in the $y,z$ plane using the formula $$ \alpha = \newcommand{\atan}{\mathrm{atan2}} \atan(u_y, u_z). $$
I am assuming here that you have the usual implementation of $\atan,$ in which $\atan(y,x)$ gives you the angle of the vector $(x,y)$ in the $x,y$ plane measured in the counterclockwise direction from the $x$ axis. But because we are doing $\atan(u_y, u_z)$ instead of $\atan(u_y, u_x),$ if you have the usual right-handed set of $x,y,z$ axes, and if you are looking onto the $y,z$ plane from the positive $x$ axis, the angle $\alpha$ is measured in a clockwise direction starting from the positive $z$ axis.
Note that the range of output of the $\atan$ function in most software implementations is either $[-\pi,\pi)$ or $(-pi,\pi].$ If $u_y > 0$ this will result in a positive angle; if $u_y < 0$ it will result in a negative angle. If $u_z > 0$ the magnitude of the angle will be less than $\pi/2,$ but if $u_x < 0$ the magnitude of the angle will be greater than $\pi/2.$ In short, no matter which way your original vector is pointing, $\atan(u_y, u_z)$ will tell you exactly how much it is rotated away from the positive $z$ direction (even if that is an obtuse angle) and in which direction, and it will identify exactly which of the four quadrants in the $y,z$ plane the vector points toward. This makes $\atan$ much more convenient than function such as $\mathrm{asin}$ or even the one-parameter $\mathrm{atan}$ function, which can only give you angles in two quadrants and therfore (in general) give you the wrong direction whenever the correct direction was in one of the other two quadrants.
If the $M_X$ parameter of your rotation function represents the conventional right-handed rotation around the $x$ axis, then $\alpha$ is exactly the rotation angle you need in order to rotate $u$ parallel to the $x,z$ plane. But you can set up a test case with an arbitrary $u$ (for example $u = (1,2,3)$) and try both rotation directions ($\alpha$ and $-\alpha$) to see which one zeros out the $y$ coordinate of the result in your software.
The result of this first rotation, viewed as a vector parallel to the $x,z$ plane, has coordinates $\left(u_x, \sqrt{u_y^2 + u_z^2}\right).$ So it makes an angle $$ \beta = \atan\left(u_x, \sqrt{u_y^2 + u_z^2} \right) $$ with the positive $z$ axis. If you are looking at the $x,z$ plane from the positive $y$ axis, the angle $\beta$ is measured counterclockwise from the positive $z$ axis, so the necessary rotation around the $y$ axis to bring the vector parallel to the $z$ axis is $-\beta$; but you can try both $-\beta$ and $\beta$ with a test vector in your software to find out which one is correct.
So the sequence of axis rotations, depending on exactly how your software interprets the axes and rotation angles, is:
$\alpha$ (or possibly $-\alpha$) around the $x$ axis;
$-\beta$ (or possibly $\beta$) around the $y$ axis.
You can then do any rotation you like (or a zero rotation) around the $z$ axis without any effect on the correctness of the result; but rotation around the $z$ axis will move both the $\startp'$ and $\endp'$ points to any direction from the $z$ axis you want. (That is, you can't change the $z$ coordinates, but you can take the $x$ and $y$ coordinates anywhere along a circle around the $z$ axis.)