Mutual information of $ \infty $-PAM I'm trying to compute using matlab the mutual information for an $ \infty $-PAM input (the limit of a very dense PAM constellation) for a range of snr and I got stuck.

I'm working with a real-valued scalar Gaussian channel with input power constrained to unity and standard Gaussian noise $ N \sim (0, 1) $. The output of the channel is given by, 
$$ y = \sqrt{ \gamma } x + n $$

I know that under the power input constraint the input is uniformly distributed on $ \left[ -\sqrt{ 3 }, \sqrt{ 3 } \right] $.

I found that in HIGH-SNR regime, there exists a gap between the Shannon capacity and the capacity achieved with $ \infty $-PAM around 0.254 bits. I know this is called the shaping loss and represents the loss of using an uniform input rather than a Gaussian input.

I know how to compute that. But, I got stuck trying to compute the capacity for a range of snr or $ \gamma $. I tried using the formula,
$$ I(x;y) = h(y) - h(n) $$
with $ f_{y}(y) = \int f_{x}(x)f_{y|x}(y|x) dx $ and I got
$$ f_{y} (y) = \frac{ 1 }{ 2 \sqrt{ 3 \gamma } } \frac{ 1 }{ 2 } \int_{ \frac{ y - \sqrt{3} }{ \sqrt{2} } }^{ \frac{ y + \sqrt{3} }{ \sqrt{2} } } \frac{ 2 }{ \sqrt{ \pi } } e^{ -t^{2} } dt \\
= \frac{ 1 }{ 2 \sqrt{ 3 \gamma } } \frac{ 1 }{ 2 } \left[ \text{erf} \left( \frac{ y + \sqrt{ 3 \gamma } }{ \sqrt{ 2 } } \right) - \text{erf} \left( \frac{ y - \sqrt{ 3 \gamma } }{ \sqrt{ 2 } } \right) \right] $$

I tried to compute $ h(y) = \int_{-\infty}^{\infty} f_{y} (y) \log_{2} f_{y} (y) $ and the integral does not converge. I tried with int$(\cdot)$ and integral$(\cdot)$ methods in matlab.

I also tried approximating the $ \text{erf}(\cdot) $ function with Taylor series expansion
$$ erf(x) = \frac{ 2 }{ \sqrt{ \pi } } \sum_{ n = 0 }^{ \infty } \frac{ (-1)^{n} x^{ 2n + 1 } }{ n!(2n + 1) } $$
but it didn't work.

How can I obtain a more manageable integral? Is my approach correct?

I appreciate any help!! Thanks.
I got an answer here, 

https://math.stackexchange.com/questions/1554659/mutual-information-for-a-continuous-uniform-distribution

The $ \text{erf}(x) $ function approximates to 1 with an error less than 0.0001 for $ x \geq 3 $. I plotted $ f_{y} (y) $ and I saw that the most of probability is within a band denpending on $ \sqrt{ 3 \gamma } $. Moreover, taking the limit $ \lim_{ y\to\infty } g(y,\gamma) \log g(y,\gamma) = 0 $. Then, I computed

$$ h(y) = \frac{1}{2} \log \left( 24 \gamma \right) - \frac{ 1 }{ 2 \sqrt{ 3\gamma} } \frac{1}{2} \int_{ - \sqrt{ 3\gamma } - a }^{ \sqrt{ 3\gamma } + a } g(y,\gamma) \log_{2}( g(y,\gamma) ) dy $$

with

$$ g (y,\gamma) = \text{erf} \left( \frac{ y + \sqrt{3\gamma} }{ \sqrt{2} } \right) - \text{erf} \left( \frac{ y - \sqrt{3\gamma} }{ \sqrt{2} } \right) $$

using $ quadgk(\cdot) $ function in matlab.