function y = myfun(x, mu, sigma)
% normal distribution function for integral computing
    y = exp(-(x-mu).^2./(2*sigma^2)) ./ (sqrt(2*pi)*sigma);
end