Approximate distribution of a complicated function of a random variable If $X$ is a random variable cdf $F(x)$ such that $F$ is invertible then we have the standard method of finding the pdf of any function of $X$, say,  $\sin(X) $  or $ X^3+1 $.However,in many situations we have very complicated composite  functions of random variables(with known distributions ,however).To take an example ,let $X_1,X_2,\dots,X_n$ be iid random variables following uniform distributions in the interval (0,1).Consider the following random varibales
$$ Y=\max_{1 \leq i
\leq n}X_i$$
$$ Z= \frac{1}{2}\left\{ 1+X_{n-1}+ \sqrt{\left(1-X_{n-1}\right)^2+4Y} \right\}$$
To the best of my knowledge, it is almost impossible to analytically find the distribution of the  $Z$. However, it is straight-forward to numerically simulate $Z$.I want to know if there are some methods of finding the approximate distributions of random variables defined like $Z$. Any references/links/hints/suggestions will be greatly appreciated. I believe the problem is of immense significance in the theory of probability and machine learning and other applied disciplines, but I honestly do not know much about it.

In this example we could just work analytically.  The joint distribution of $(X_{n-1},Y)$ can be written down explicitly ($Y=\max X_{n-1},Y^*$,where $Y^*$ is the maximum of the other $n-1$ Xs). We can also write down the partial derivatives of $Z$ with respect to $X$ and $Y$, and work out how many $(X,Y)$ pairs map to each $Z$. It would be sufficiently annoying that I'm not going to do it, but it doesn't seem infeasible.
Or, if $n$ is large, we can note that $X_{n-1}$ will be approximately independent of $Y$, and find the analytic solution for this simplified problem.
A solution that's more generally applicable would be to use simulation.  Simulate the $X_i$, or just simulate $(X_{n-1}, Y)$, compute $Z$, and use a kernel density estimate or a histogram or whatever to represent the pdf of $Z$ to as high accuracy as you want.
Here's some completely unoptimised code that simulates a million $Z$s in about ten seconds
sim_z<-function(){
  n<-42
  xi<-runif(n)
  Y<-max(xi)
( 1/2)*(1+xi[n-1]+sqrt((1-xi[n-1])^2+4*Y))
}


zs<-replicate(1e6,sim_z())


You could then try some partial analytic arguments: is the maximum exactly at 2 (yes)?  Does the shape of the lower tail have a simple relationship to the shape of the lower tail of $Y$? Can we work out exactly where the peak is (or approximately, for large $n$)? And so on.
