How to calculate the integral in normal distribution?

There are several methods to approach this, but I am going to use one that meets your requirement (clarified in a comment) that one must forego the use of computational engines like Mathematica, instead opting for a calculator. Additionally, I feel that using a table of normal distribution values is cheating, so I will be foregoing their use as well.

First, we need the equation for $\mathcal{N}(0,25)$, which, by definition, is: \begin{align*} f(x) &= \mathcal{N}(\mu,\sigma ^2)\\ &= \mathcal{N}(0,25)\\ &= \frac{1}{\sigma\sqrt{2\pi}}\,e^{ -\frac{(x-\mu)^2}{2\sigma^2} }\\ &= \frac{1}{5 \sqrt{2 \pi }}\,e^{-\frac{x^2}{50}} \end{align*} Now, we simply need to integrate this from $-x$ to $x$, set it equal to $.90$, and solve for $x$ (our answer): $$F(x) = \frac{1}{5 \sqrt{2 \pi }}\int_{-x}^x e^{-\frac{x^2}{50}} \,\mathrm{d}x=0.9$$ However, we run into problems when we realize that this isn't a simple integral to take! Since we already know the answer, we can actually exploit this to our advantage by finding an elementary function that estimates $F(x)$ with an adequately small margin of error after the integral. One simple function you can use to estimate $F(x)$ is a Taylor series.

We first need to find a Taylor series for $f(x)$ using the formula for a Taylor series: $$\sum\limits_{n=0}^{\infty } \frac {f^{(n)}(a)}{n!} \, (x-a)^{n}$$ One can easily recognize the pattern for our function when $a=0$ (the center of our bell curve) to generate this series: $$f(x) = \frac{1}{5 \sqrt{2 \pi}}\sum\limits_{k=0}^{\infty} \frac{(-1)^k x^{2k}}{50^k k!}$$ We can now trivially take the integral of this series where we otherwise wouldn't have been able to: $$F(x) = \frac{1}{5 \sqrt{2 \pi}}\sum\limits_{k=0}^{\infty} \frac{(-1)^k x^{2 k+1}}{50^k k! (2 k+1)}$$ Since we can't use our pencil and paper to evaluate the series all the way to infinity, we need to figure out how many terms we need to go out to get an acceptably accurate answer. This can be accomplished by figuring out the error caused by not going out to infinity before we actually solve the equation. By plugging in our already-known value of $8.225$ for $x$ in the $k$th term of the series and evaluating just the series part from $k$ to $\infty$, we get the error that the series will have from $0$ to $k-1$. Since we, again, can't go all the way out to infinity, we can get a slight, but adequate, underestimate of the error by simply evaluating the series from $k$ to $k$ since each subsequent term in the series is substantially smaller.

I started by plugging in $8.225$ for $x$ when $k=7$ and got this (do this on your calculator): $$\sum\limits_{k=7}^{7 } \frac{(-1)^k (8.225)^{2 k+1}}{50^k k! (2 k+1)} = -\frac{(8.225)^{15}}{59062500000000000} \approx -0.000903081$$ Since our integrated series is the equivalent of $F$ in $F(x) - F(-x) = \int_{-x}^x f(x)\,\mathrm{d}x$ where $f(x)$ is the equation for $\mathcal{N}(0,25)$, we need to multiply our answer by $2$ to get a rough estimate of the error our final answer will have: $-0.000903081 \times 2 = -0.001806162$.

Since we are looking for a three decimal accuracy, we should proceed and try the series when $k=8$: $$\sum\limits_{k=8}^{8 } \frac{(-1)^k (8.225)^{2 k+1}}{50^k k! (2 k+1)} = \frac{(8.225)^{17}}{26775000000000000000} \approx 0.000134766$$ Finally, $0.000134766 \times 2 = 0.000269532$. This is of adequate accuracy for our final calculation; since the error calculated when $k=8$ is acceptable, we will evaluate the terms of the series from $k=0$ to $7$ ($k$ used in error calculation minus $1$): $$F(x) \approx \frac{1}{5 \sqrt{2 \pi}}\sum\limits_{k=0}^{7} \frac{(-1)^k x^{2 k+1}}{50^k k! (2 k+1)}$$ This comes out to be: $$\frac{x-\frac{x^3}{150}+\frac{x^5}{25000}-\frac{x^7}{5250000}+\frac{x^9}{1350000000}-\frac{x^{11}}{412500000000}+\frac{x^{13}}{146250000000000}-\frac{x^{15}}{59062500000000000}}{5 \sqrt{2 \pi }}$$ Since this is our estimate of $F(x)$, and this is an odd function (a function is odd if $-f(x) = f(-x)$), we simply need to multiply this function by $2$ to get $F(x) - F(-x)$ (what we need to set to $0.90$ and solve): $$\frac{1}{5} \sqrt{\frac{2}{\pi}} (x-\frac{x^3}{150}+\frac{x^5}{25000}-\frac{x^7}{5250000}+\frac{x^9}{1350000000}-\frac{x^{11}}{412500000000}+\frac{x^{13}}{146250000000000}-\frac{x^{15}}{59062500000000000})$$

Now, by setting this equal to $.90$, rearranging the equation as a polynomial, and using a method of our choice to solve polynomials on a calculator (like Newton's method to converge on the answer), we find that the relevant solution is: $$x \approx 8.22473 \approx 8.225$$


There are several ways to compute the cumulative normal distribution.


Simple Series Integration

First of all, we can start with $$ e^{-x^2/2}=1-\frac{x^2}{2^1\cdot1!}+\frac{x^4}{2^2\cdot2!}-\frac{x^6}{2^3\cdot3!}+\dots $$ and integrate to get $$ \begin{align} \frac1{\sqrt{2\pi}}\int_0^xe^{-t^2/2}\,\mathrm{d}t &=\frac1{\sqrt{2\pi}}\left(x-\frac{x^3}{3\cdot2^1\cdot1!}+\frac{x^5}{5\cdot2^2\cdot2!}-\frac{x^7}{7\cdot2^3\cdot3!}+\dots\right)\\ &=\frac1{\sqrt{2\pi}}\sum_{k=0}^\infty(-1)^k\frac{x^{2k+1}}{(2k+1)2^kk!} \end{align} $$


Unilateral Power Series

To get a series with no sign changes, write $$ \Omega(x)=e^{x^2/2}\int_0^xe^{-t^2/2}\,\mathrm{d}t $$ and note that $$ \Omega'(x)=1+x\Omega(x) $$ which leads to the following recursion for the coefficients $$ a_{n+2}=\frac{a_n}{n+2} $$ Because $\Omega(0)=0$ and $\Omega'(0)=1$, we get $$ \begin{align} \frac1{\sqrt{2\pi}}\int_0^xe^{-t^2/2}\,\mathrm{d}t &=\frac1{\sqrt{2\pi}}e^{-x^2/2}\Omega(x)\\ &=\frac1{\sqrt{2\pi}}e^{-x^2/2}\sum_{k=0}^\infty \frac{x^{2k+1}}{(2k+1)!!} \end{align} $$


Asymptotic Expansion

The terms of the convergent series above can get quite big as $x$ gets big, so let's compute an asymptotic expansion for the tail.

Consider $$ \begin{align} \int_x^\infty e^{-t^2/2}\,\mathrm{d}t &=\int_0^\infty e^{-(t+x)^2/2}\,\mathrm{d}t\\ &=e^{-x^2/2}\int_0^\infty e^{-xt-t^2/2}\,\mathrm{d}t\\ &=\frac{e^{-x^2/2}}{x}\int_0^\infty e^{-u-u^2/(2x^2)}\,\mathrm{d}u\\ &=\frac{e^{-x^2/2}}{x}\int_0^\infty e^{-v}u'\,\mathrm{d}v\\ \end{align} $$ where $v=u+\dfrac{u^2}{2x^2}$; that is, $u=x^2\left(\sqrt{1+2v/x^2}-1\right)$. Using the binomial theorem, to expand $\sqrt{1+2v/x^2}$, we get $$ \frac1{\sqrt{2\pi}}\int_x^\infty e^{-t^2/2}\,\mathrm{d}t \sim\frac1{\sqrt{2\pi}}e^{-x^2/2}\sum_{k=0}^\infty\frac{(-1)^k(2k-1)!!}{x^{2k+1}} $$ where $(-1)!!=1$.

Note that this series is not convergent, but asymptotic.