It took a while, but I have an answer.
This answer describes the transformation of a hyperbolic paraboloid into canonical form. The procedure is similar for other quadric shapes.
An affine transform $\mathbf{T}$ is of the form:
\begin{equation}
\mathbf{T} = \begin{bmatrix}
T_{11} & T_{12} & T_{13} & T_{14} \\
T_{21} & T_{22} & T_{23} & T_{24} \\
T_{31} & T_{32} & T_{33} & T_{34} \\
0 & 0 & 0 & 1
\end{bmatrix}
\end{equation}
Every quadric may be transformed via affine transformations into either a diagonal matrix or a near diagonal matrix, i.e.
\begin{equation}
\mathbf{X}^T \mathbf{T}^T \mathbf{A} \mathbf{T} \mathbf{X} = \mathbf{X}^T \mathbf{D} \mathbf{X}
\end{equation}
where $\mathbf{D}$ is a (near-)diagonal matrix. A near-diagonal matrix is symmetric and is of the form:
\begin{equation}
\mathbf{D} = \begin{bmatrix}
D_{11} & 0 & 0 & 0 \\
0 & D_{22} & 0 & 0 \\
0 & 0 & 0 & D_{43} \\
0 & 0 & D_{43} & 0
\end{bmatrix}
\end{equation}
The canonical form of a hyperbolic paraboloid listed in Xu (Table 1, pp518) is not strictly in near-diagonal form:
\begin{align}
0 &= z - xy \\
&= \mathbf{X}^T \mathbf{D'} \mathbf{X} \\
&= \mathbf{X}^T \begin{bmatrix}
0 & -0.5 & 0 & 0 \\
-0.5 & 0 & 0 & 0 \\
0 & 0 & 0 & 0.5 \\
0 & 0 & 0.5 & 0
\end{bmatrix} \mathbf{X}
\end{align}
but may be re-written in near-diagonal form as:
\begin{align}
0 &= -x^2 + y^2 + z \\
&= \mathbf{X}^T \mathbf{D} \mathbf{X} \\
&= \mathbf{X}^T \begin{bmatrix}
-1 & 0 & 0 & 0 \\
0 & 1 & 0 & 0 \\
0 & 0 & 0 & 0.5 \\
0 & 0 & 0.5 & 0
\end{bmatrix} \mathbf{X}
\end{align}
The Xu form $\mathbf{D'}$ and the near-diagonal form $\mathbf{D}$ are related via an affine transform:
\begin{align}
\mathbf{D'} = \mathbf{M}^T \mathbf{D} \mathbf{M}
\end{align}
where:
\begin{equation}
\mathbf{M} = \begin{bmatrix}
0.5 & -0.5 & 0 & 0 \\
0.5 & 0.5 & 0 & 0 \\
0 & 0 & 1 & 0 \\
0 & 0 & 0 & 1 \\
\end{bmatrix}
\end{equation}
Since the Xu's canonical form can be found from the near-diagonal form, our immediate task is to find the transformation $\mathbf{T}$ that casts a given hyperbolic paraboloid $\mathbf{A}$ into near-diagonal form. Although Dupont recommends determining the transform via Gaussian Reductions to avoid radicals, perhaps the easiest first step is to remove the rotation via an Eigendecomposition as suggested by Levin.
Consider the principle submatrix $\mathbf{A}_u$. Since it is symmetric, of rank 2 and has known signs for the eigenvalues, an eigendecomposition may be found such that
\begin{align}
\mathbf{A}_u = {\mathbf{R}_u}^T \mathbf{\Lambda} \mathbf{R}_u
\end{align}
where $\mathbf{R} \in SO(3)$ and $\mathbf{\Lambda}$ is in the form:
\begin{equation}
\mathbf{\Lambda} = \begin{bmatrix}
\lambda_+ & 0 & 0 \\
0 & \lambda_- & 0 \\
0 & 0 & 0 \\
\end{bmatrix}
\end{equation}
where $\lambda_+ > 0$ and $\lambda_- < 0$. Therefore, the first affine transform may be built as:
\begin{align}
\mathbf{T}_1 = \begin{bmatrix}
\mathbf{R}_u & \mathbf{0} \\
\mathbf{0} & 1
\end{bmatrix}
\end{align}
The second transform forces the eigenvalues to unit magnitude and it can be verified that:
\begin{align}
\mathbf{T}_2 = \begin{bmatrix}
\frac{1}{\sqrt{\lambda_+}} & 0 & 0 & 0 \\
0 & \frac{1}{\sqrt{\|\lambda_-\|}} & 0 & 0 \\
0 & 0 & 1 & 0 \\
0 & 0 & 0 & 1 \\
\end{bmatrix}
\end{align}
Let $\mathbf{A}_2 = {\mathbf{T}_2}^T {\mathbf{T}_1}^T \mathbf{A} \mathbf{T}_1 \mathbf{T}_2$, which will be of the form:
\begin{equation}
\mathbf{A}_2 = \begin{bmatrix}
1 & 0 & 0 & a \\
0 & -1 & 0 & b \\
0 & 0 & 0 & c \\
a & b & c & d
\end{bmatrix}
\end{equation}
It can be verified that the transform required to cast into near-diagonal form is:
\begin{equation}
\mathbf{T}_3 = \begin{bmatrix}
1 & 0 & 0 & 0 \\
0 & 1 & 0 & 0 \\
\frac{-a}{c} & \frac{-a}{b} & \frac{-1}{2c} & \frac{-d}{2c} \\
0 & 0 & 0 & 1
\end{bmatrix}
\end{equation}
Therefore, the complete transformation required to cast the hyperbolic paraboloid into Xu's canonical form is:
\begin{align}
\mathbf{T} = \mathbf{T}_1 \mathbf{T}_2 \mathbf{T}_3 \mathbf{M}
\end{align}
Your prescription sounds right to me.
Suppose $R=LL'$ is the desired covariance matrix with its Cholesky decomposition. Let e be a random (column) vector distributed as a zero-mean multivariate Gaussian with unit variance and no correlations. That is, $E[e]=0$ and $E[ee'] = I$. Then $L$ and $e$ are the results of your steps 1 and 2.
The matrix-vector product $x=Le$ is the result of your step (3). I'll assert that $x$ is also a multivariate Gaussian, since proof of that is not the issue here. We can see that the expected covariance of $x$ is
$$\text{Covar}(x) = E[xx']-E[x]E[x'] = E[ L e e' L'] - 0 = L E[ee'] L' = LIL' = LL' = R.$$
So $x$ should follow a multivariate Gaussian distribution with the correct covariance matrix, $R$. All I can think of is that if you want to error ellipse to be distributed about a point other than the origin, then (step (4)) you need to add that vector to $x$.
Oh, I know one potential problem: Perhaps you used the upper, rather than the lower Cholesky factor? If you let $x=L'e$, then your random vectors would have an expected covariance of $L'L$, which does not equal $R$ unless $R$ is diagonal (i.e., this only works if there's no correlation, which is surely not the case you have in mind).
Best Answer
TL;DR: Just add zeros.
The general quadric equation can be broken down into the sum of a quadratic form, a linear term and a constant: $$\mathbf x^T\mathtt A\mathbf x+2\mathbf b^T\mathbf x+c=0.$$ The coefficients can be packaged up into a single homogeneous matrix $\mathtt Q$ by setting $\hat{\mathbf x}=(\mathbf x^T; 1)^T$ so that the equation becomes $$\hat{\mathbf x}^T\mathtt Q\hat{\mathbf x} = \begin{bmatrix}\mathbf x&1\end{bmatrix} \begin{bmatrix}\mathtt A & \mathbf b \\ \mathbf b^T & c\end{bmatrix} \begin{bmatrix}\mathbf x\\1\end{bmatrix} = 0.$$ (This equation is homogeneous in the components of $\hat{\mathbf x}$, so we can relax the restriction that the last component be $1$.)
For your covariance matrix, the ellipsoid equation is $\mathbf x^T\mathtt M\mathbf x=1$, therefore $\mathbf b=0$ and $c=-1$, i.e., $$\mathtt Q = \begin{bmatrix}\mathtt M & \mathbf 0 \\ \mathbf 0^T & -1\end{bmatrix}.$$ The projection matrix $\mathtt P$ for orthogonal projection onto the $x$-$y$ plane simply deletes the third component of a vector, so the matrix $\mathtt C$ that represents the projected outline of the ellipsoid will also have this form, but with a $2\times2$ matrix $\mathtt M'$. Computing the matrix of the outline thus comes down to inverting the upper-left $2\times2$ submatrix of $\mathtt M^{-1}$ and the Cartesian equation of the resulting ellipse is $(x,y)\,\mathtt M'(x,y)^T=1$. If you’re so inclined, you could even develop a closed-form formula for the coefficients of this equation from $\mathtt M^{-1}=\operatorname{cof}(M)^T/\det(M)$. The resulting expressions turn out to be quite simple.