The problem is, the distribution of $S_n$ is not normal, nor closely approximated by a normal random variable, since $3$ rolls is a small number. Note that if the first two dice give a total of at most $3$ (which can happen in $3$ ways), then no matter what the third die shows, the total will be at most $9;$ this gives $3\cdot 6=18$ ways to get a total of at most $9$ if the first two dice total at most $3$. If the first two dice give a total of $4$ (which can happen in $3$ ways), then we must get at most $5$ on the third die; this gives $3\cdot5=15$ ways to get a total of at most $9$ if the first two dice give a total of $4.$ Similarly, there are (respectively) $4\cdot4=16,$ $5\cdot3=15,$ $6\cdot2=12,$ and $5\cdot1=5$ ways to get a total of at most $9$ if the first two dice give a total of (respectively) $5,6,7,$ and $8.$ There are no other totals on the first two dice that will allow us to get a total of at most $9$ on all three dice. Hence, adding them up, there are $81$ ways to get a total of at most $9$ on the three dice. Since there are $216$ ways to roll, it follows that the probability of a total of more than $9$ is $\frac58=0.625.$ (Why?)
I leave it to you to fill in the details, but let me know if you have any questions.
We are trying to use Monte Carlo Simulation to find:
$$\displaystyle \int_0^\infty \dfrac 34 x^4e^{-x^{3/4}}\,dx$$
As you have discovered, the infinite limits are problematic. We can do several things to overcome this
Approach 1: A naive approach is to plot the function given that it has a decaying exponential which tends to zero very quickly. If we plot this function, we see:
So, lets try to approximate it as:
$$\displaystyle \int_0^{50} \dfrac 34 x^4e^{-x^{3/4}}\,dx$$
Using your approach (recall, now your code needs to change to generate random numbers $(a, b)$ and with the function of interest) with $n = 10000$, we get:
$$I = 387.256$$
If we calculate the exact integral we get:
$$I = \Gamma \left(\frac{20}{3}\right) \approx 389.035$$
Approach 2: In this approach we will make a change in variables so as to map an infinite range of integration into a finite one.
A good example that works for functions that decrease towards infinity faster than $\frac{1}{x^2}$ is to use:
$$\displaystyle \int_a^b f(x)~dx = \int_{\frac{1}{a}}^{\frac{1}{b}} \dfrac{1}{t^2}f\left(\dfrac{1}{t}\right)~dt$$
With this altered approach, we can now do integrals when one of the limits is infinite
and the other limit is of the same sign. If you need to integrate say something that has $a = 0$, what you do is break the integral into two parts and use Monte Carlo on each one.
If we have infinite limits for both limits, we could split the integral into three or more parts. So, for this problem, we have to split it into two pieces as:
$$\displaystyle \int_0^\infty \dfrac 34 x^4e^{-x^{3/4}}\,dx = \int_0^1 \dfrac 34 x^4e^{-x^{3/4}}\,dx + \int_0^1 \dfrac{1}{t^2} \left(\dfrac 34\right) \left(\dfrac 1t \right)^4 e^{{-(\frac 1t)}^{\frac 34}}~dt$$
Using the Monte Carlo for the first part, we get $I_1 = 0.0647226$ and for the second part, we get $I_2 = 387.551$, which yields:
$$ \displaystyle \int_0^\infty \dfrac 34 x^4e^{-x^{3/4}}\,dx = 387.616$$
Compare that to our previous and the exact result.
Approach 3: In this approach we will make a change in variables so as to map an infinite range of integration into a finite one. This time, we will transform variables using $u = \frac{1}{1+x}$, yielding:
$$\displaystyle \int_0^\infty \dfrac 34 x^4e^{-x^{3/4}}\,dx = \int_0^1 \frac{3 e^{-\left(\frac{1-u}{u}\right)^{3/4}} (1-u)^4}{4 u^6} ~ du$$
Running the Monte Carlo simulation with $n = 10000$ yields:
$$I = 390.18704821709946$$
Best Answer
Electro82 was definitely on the right track - just a couple of mods here and this works: