The problem of determining distributions of the product of random variables is one of the most important problems. However, most studies only focus on independence structures on some common distributions of the functions of random variables or stochastic dependence through multivariate normal joint distributions or unknown joint distribution. To bridge the gap in the literature, in this paper we first derive the general formulas to determine the density and distribution for the product of two or more random variables via copulas to capture the dependence structures among the variables. We then propose an approach combining Monte Carlo algorithm, graphical approach, and numerical analysis to efficiently estimate both density and distribution. Thereafter, we illustrate our approach by examining the shapes and behaviors of both density and distribution of the product of two log-normal random variables on several different copulas, including Gaussian, Student-t, Clayton, Gumbel, Frank, and Joe Copulas, and estimates some common measures including mean, median, standard deviation, skewness, and kurtosis. We find that different types of copulas strongly affected behavior of distributions differently. For example, we find that the product is strongly affected in median, variance, skewness, and kurtosis when using copulas from the elliptical family but not strongly affected when using copulas from the Archimedean family. We have discussed the behaviours of all copulas with the same Kendall coefficient and drawn conclusions on our findings. Our results are the foundations of any further study that relies on the density and cumulative probability functions of product of n random variables and the theory we developed in this paper is useful to both academics, practitioners, and policy makers.