Initially, there is a random variable $X$ (non-negative) with distribution function: $$P(x) = \lambda e^{-\lambda x}$$ Now we randomly generate $X$ and, then, form sets of $N$ values of this variable: $(x_1,x_2,\ldots,x_N)$
And for each set, we define, for example, $$S = \frac{1}{N}\sum_{k=1}^N x_k^2$$ $S$ is then a new random variable. My question is: How to find the pdf of this new random variable?
Comments; Your Question is a possible duplicate of this Q&A, where a method is shown to get the answer, but not the answer itself.
You mention simulation: Here are results from a simulation in R of 100,000 realizations of $X^2,$ where $X \sim \mathsf{Exp}(\text{rate}\, \lambda = 1/3)$ and $n = 5.$ Your reason for simulation is unclear; perhaps you want to see if theoretical and simulated answers agree.
Of course, if $n$ is large enough, this sum of independent random variables will be approximately normal. With $m = 100000$ the mean and SD should be accurate to a couple of significant digits.