Suppose that $X_1 ,\ldots,X_n,Y_1,\ldots,Y_n$ are all independent normal random variables with different means and variances. What is the PDF of the following random variable?
$$X_1Y_1+\cdots+X_nY_n$$
or is there any way that I can find an approximation for its PDF?
On
You want to know the distribution of
$$ Z = \sum_{i=1}^n X_iY_i, $$
where $X_i$ and $Y_i$ are all normal and independent but not identically distributed. Let's calculate the mean and variance of $Z$:
$$ \mathrm E[Z] = \sum_{i=1}^n \mathrm E[X_iY_i] = \sum_{i=1}^n \mathrm E[X_i]\mathrm E[Y_i]. $$
You know the individual means of the normals, so the above sum can be calculated. Now the variance, knowing that covariance terms drop out because of uncorrelation due to independence
$$ \mathrm{Var}[Z] = \sum_{i=1}^n \mathrm{Var}[X_iY_i] = \sum_{i=1}^n \mathrm E[Y_i]^2\mathrm{Var}[X_i]+\mathrm E[X_i]^2\mathrm{Var}[Y_i]+\mathrm{Var}[X_i]\mathrm{Var}[Y_i] $$
Again, this is all in terms of quantities you know, the means and variances of the underlying distributions.
So, IF the CLT applies, your answer is that $Z\sim\mathcal{N}(\mathrm{E}[Z],\mathrm{Var}[Z])$. Rigorously proving that this is true for your particular case and your particular distributions is up to you, but my claim is that it is true. I pastebinned some numerical calculations to support my claim.
The product $Z_i=X_iY_i$ obeys the product normal distribution, which has a characteristic function given in: Characteristic function of product of normal random variables. As a result the characteristic function of your random variable is given by a product of the similar characteristic function (after matching to the various means and standard deviations). You can solve the question by approximating the Fourier transform of the resulting characteristic function.