My problem is the following. I am given a PDF, say $p(x) = x \cdot \mathrm{e}^{-x}$ for $x > 0$. I want to generate $n$, say $n = 1000$, random variables given this density, such that I can calculate the probability that the $Y = X_{1} + \cdots + X_{n} > \gamma$ where $\gamma$ is large. How do I sample these random variables given the PDF. Do I need to calculate the mean and variance of $Y$ and then draw these $X$'s from a normal distribution with $\mu_{Y}$ and $\sigma_{Y}^{2}$. I would like to implement this into Python.
Btw. I have tried to integrate $p(x)$ to get its CDF and then use the inverse method, but unfortunately it does not work.
In general you would have to numerically invert the cdf as described in comments, but the specific pdf of $xe^{-x}$ gives the gamma distribution with shape parameter $k=\alpha=2$. This distribution has the nice property that the sum of same-scale (including iid) random variables is also gamma-distributed: $$\sum_i\Gamma(k_i,\theta)\sim\Gamma\left(\sum_ik_i,\theta\right)$$ So your $Y$, being the sum of $n$ iid $\Gamma(2,1)$ variables, is $\Gamma(2n,1)$ distributed. Then you can calculate $P(Y>\gamma)$ from there.