Your normal distribution is one that is too narrowly defined around the mean of 350. In fact, it is effectively zero towithin tolerances outside of 350 +/- 60, as that would be 6 sigma in each direction. So what happens is for most limits of integration, quadgk sees something that seems to be zero everywhere it looks. Therefore the integral must be zero.
That is, over the domain [0,1000], that kernel can look to be an essential Dirac delta, centered at 350.
This is a common mistake made by people. When I see someone claim that a numerical integration routine is incorrect, I first would look to see if they have a narrow normal distribution in the kernel.
So, only if I look in a narrow region do I see anything that is even slightly different from 0. What is the integral of a function that is 0 over any domain? Zero. You always need to be careful when numerically integrating a function that looks like a delta function. And normal distributions tend to look exactly like that when viewed over a too wide interval.
This is the code you wrote:
Ewr=@(mu) quadgk(@(x) f1(x,mu),0,(mu-1/wo)*T)
So, f1 is not a function of mu, which enters into the problem only as the upper limit of integration. Essentially, you are varying B from 0:1000, which then appears as mu.
But now look closely at what you do. I'll just use a limited set of values here.