How do I find the probability involving the normal distribution function for $P(X \le 16.50)$?

Here $\sigma=1.01$, $\mu=19.47$, s.t., we have

$P(X\leq 16.50)$

$= \frac{1}{\sigma\sqrt{2\pi}}\int\limits_{-\infty}^{16.50}e^{-\frac{{(x-\mu)^2}}{2\sigma^2}}dx$

$= \frac{1}{\sigma\sqrt{2\pi}}\int\limits_{-\infty}^{(16.50-19.47)/1.01}e^{-\frac{z^2}{2}}\sigma dz$, substitute $z=\frac{x-\mu}{\sigma}$, s.t., $dx=\sigma dz$

$= \frac{1}{\sqrt{2\pi}}\int\limits_{-\infty}^{-2.940594}e^{-\frac{z^2}{2}}dz$

$= \frac{1}{\sqrt{2\pi}}\int\limits_{-\infty}^{0}e^{-\frac{z^2}{2}}dz-\frac{1}{\sqrt{2\pi}}\int\limits_{-2.940594}^{0}e^{-\frac{z^2}{2}}dz$

$= \frac{1}{2}-\frac{1}{\sqrt{2\pi}}\int\limits_{-2.940594}^{0}e^{-\frac{z^2}{2}}dz$, being standard normal pdf symmetric around $z=0$

$= \frac{1}{2}-\frac{1}{\sqrt{2\pi}}\int\limits_{-2.940594}^{0}\left(1-\frac{z^2}{2}+\frac{z^4}{8}-\frac{z^6}{48}+\frac{z^8}{384}-\frac{z^{10}}{3840}+\ldots\right)dz$, by Maclurin

$= \frac{1}{2}-\frac{1}{\sqrt{2\pi}}\left[z-\frac{z^3}{6}+\frac{z^5}{40}-\frac{z^7}{336}+\frac{z^9}{3456}-\frac{z^{11}}{42240}+\ldots\right]_{-2.940594}^{0}$, compute this sum for at least $20$ terms (the more terms you consider the better precision you get)

$=0.5-\frac{1.249209}{\sqrt{2\pi}}=0.5-0.4983623 = 0.0016377$

You can use the following R code to compute the sum of first $20$ terms of the series:

z <- -2.940594
s <- z
for (i in 1:20) {
   s <- s + (-1)^i*z^(2*i+1)/((2*i+1)*2^i*factorial(i))
}
s
# -1.249209

Let's verify the result by numeric integration:

integrate(function(z) exp(-z^2/2), lower=-Inf, upper=-2.940594)$value / sqrt(2*pi)
#[1] 0.001637919

and by computing the normal CDF with pnorm():

pnorm(16.50, mean = 19.47, sd = 1.01)
#[1] 0.001637918