On solving standard normal deviation equation numerically Let $\mu$, $\sigma$, and $a$ be given real numbers and $\epsilon>0$ given. What kind of methods there are to solve the equation $$\dfrac{1}{\sqrt{2\pi}\sigma}\int_{-\infty}^x e^{-\dfrac{(t-\mu)^2}{2\sigma^2}}dt=a$$ numerically such that error is at most $\epsilon$? In particular, is it easier to solve the equation by converting it to the standard normal distribution rather than approximating exp to its series expansion and solve polynomial equation numerically? At least the problem seems hard if $\sigma \approx 0$. 
 A: Based on the Gauss-Hermite polynomial chaos expansion of a general normal random variable $X \sim \mathcal{N}(\mu,\sigma^2)$ about a standard normal (Gaussian) random variable $Y \sim \mathcal{N}(0,1)$,
$$X = \mu + \sigma Y.$$
Therefore, the most accurate way to compute the integral is to compute the appropriate integral for the standard normal distribution; any exigent numerical error will be due to a single multiplication and addition, which is bound by the machine epsilon.

As far as why I invoke the Gauss-Hermite chaos, if you use a simple change of variable argument, then you have a seeming cart-before-the-horse scenario, where the change of variable argument leads you to the same question in the OP, so there appears to be little assurance that the computation would be bounded in such a way, particularly when $\sigma \approx 0$.
However, it can be rigorously shown that the coefficients in the Gauss-Hermite chaos expansion $X = \sum_{i=0}^{\infty} k_i \Phi_i(Y)$, where $\Phi_i$ is the $i^{\textrm{th}}$ Hermite polynomial about $Y$, are all zero when $i \ge 2$, and that $k_0 = \mu$ and $k_1 = \sigma$.
This avoids the ambiguity, and provides a constructive means by which to compute the integral -- at least in-so-much as one can accurately establish the integration bounds.
