The CLT can be stated in different ways. By generating samples of different sizes and calculating the sample mean for each sample, we established that the sample mean is normally distributed with a mean equal to the population mean and a standard deviation (uncertainty) = population standard deviation divided by square root of N (that is the sample size).
I don’t see why Monte-Carlo can’t be use for demonstrating CLT. The proof might not be vigorous, but it serves the purpose, from a data science point of view. Data science is not mathematics, and mathematics is not data science.
Thanks for your criticism and feedback!