The co-variance matrix of any random vector $Y$ is given as $\mathbb{E} \left(YY^T \right)$, where $Y$ is a random column vector of size $n \times 1$. Now take a random vector, $X$, consisting of uncorrelated random variables with each random variable, $X_i$, having zero mean and unit variance $1$. Since $X_i$'s are uncorrelated random variables with zero mean and unit variance, we have $\mathbb{E} \left(X_iX_j\right) = \delta_{ij}$. Hence, $$\mathbb{E} \left( X X^T \right) = I$$ To generate a random vector with a given covariance matrix $Q$, look at the Cholesky decomposition of $Q$ i.e. $Q = LL^T$. Note that it is possible to obtain a Cholesky decomposition of $Q$ since by definition the co-variance matrix $Q$ is symmetric and positive definite.
Now look at the random vector $Z = LX$. We have $$\mathbb{E} \left(ZZ^T\right) = \mathbb{E} \left((LX)(LX)^T \right) = \underbrace{\mathbb{E} \left(LX X^T L^T\right) = L \mathbb{E} \left(XX^T \right) L^T}_{\text{ Since expectation is a linear operator}} = LIL^T = LL^T = Q$$ Hence, the random vector $Z$ has the desired co-variance matrix, $Q$.
Consider the covariance matrix of $(X,Y,Z)^T$, $C=\mathbb E[[X,Y,Z]^T[X,Y,Z]]$. (since the means are $0$)
Compute the individual variance and standard deviation of $X,Y,Z$, call it $\sigma_X,\sigma_Y,\sigma_Z$.
Define the covariance of normalized variables $\frac{X}{\sigma_X},\frac{Z}{\sigma_Z},\frac{Z}{\sigma_Z}$, we call it $\bar C$.
For this normalized covariance matrix, it shall has this structure
$$\bar C=
\begin{bmatrix}
1 & 0.2 & \rho\\
0.2 & 1 & 0.2\\
\rho & 0.2 & 1\\
\end{bmatrix}
$$
Note scale and shifting of random variables do not change the pearson correlation among them. so here $\rho=\rho_{XZ}$
To make $\bar C$ a valid covariance matrix, it needs to be positive (semi-)definite. Using a generalized version of Sylvester criterion for pos-semi-definite matrices, we need to show all principal minors are non-negative (in contrast to leading ones are positive).
This criterion comes down to
$$
\det(\bar C)\geq 0\\
\det(\begin{bmatrix}
1 & \rho\\
\rho & 1\\
\end{bmatrix})\geq0
$$
Which translates to
$$
1-\rho^2\geq 0\\
\frac{576}{25^2}-(\rho-\frac 1{25})^2\geq 0
$$
Then you get the range of the correlation between $X,Z$. $[-\frac{23}{25},1]$
$$
-\frac{23}{25}\leq \rho\leq 1
$$
For the last part, notice how linear transform affects the covariance of normal variables. Generally, for multivariate Gaussian
$U\sim \mathcal N(0,C)$ the transform it with an invertible transform $V=AU$, $V\sim \mathcal N(0,ACA^T)$.
Thus if you have three independent standard normal variable $[N_1,N_2,N_3]^T\sim \mathcal N(0,I)$. Then if we have a Cholesky decomposition of $C=LL^T$ we can define the random variables $L[N_1,N_2,N_3]^T\sim \mathcal N(0,LL^T)=\mathcal N(0,C)$
Thus we find the transform that maps independent normal variables to $X,Y,Z$
$$
[X,Y,Z]^T=L[N_1,N_2,N_3]^T
$$
Prussing, John E., The principal minor test for semidefinite matrices, J. Guid. Control Dyn. 9, 121-122 (1986). ZBL0599.15010.
Best Answer
From a matrix algebra point of view the answer is fairly simple. Assume your covariance matrix is $\Sigma$ and let
$$ D =\sqrt{ \text{diag}\left( {\Sigma} \right)} $$
then the correlation matrix is given by $$ \varrho = D^{-1}\Sigma D^{-1} $$
Edit: fixed to include square root