import numpy as np


def f(x, m_u, s_igma):
    denuminator = s_igma*np.sqrt(2*np.pi)
    numinator = np.exp(-0.5*((x - m_u)/sigma)**2)
    return numinator/denuminator

mu = 0.0
sigma = 1.0
num_std = 2 #number of standard dev, sigma, off the mean


start = mu - num_std*sigma
stop = mu + num_std*sigma
n = 1000 #number of subintervals

dx = (stop - start)/n #length of each subinterval

x=np.linspace(start, stop, n+1)

area=0.5*(f(start, mu, sigma) + f(stop, mu, sigma))*dx

for i in range(1, n, 1):      
    area = area + f(x[i], mu, sigma)*dx

print(round(area,3), '% of the values lie within',
      num_std, 'standard deviations.')