The problem is that you are integrating an oscillatory function over an infinite interval. The MATLAB website doesn't give specifics on the algorithm behind their
integral function (it just says 'globally adaptive'). The older
quad function used to be adaptive Simpson, so I can assume that
integral is the same. What probably is happening is that the code transforms the interval $[0,\infty]$ to $[0,1]$ or $[-1,1]$ or so and hence it is squeezing all the oscillations in this interval. As a result it is adding a lot of positive contributions and subtracting a lot of negative contributions. Hence a totally wrong answer.
For oscillatory integrals like these a couple of tricks exist. First of all, if the function is rapidly decaying, you could truncate the interval (and as such the classical routines will not transform the interval). If you have ideas on the zero crossings, you could loop over all these and integrate between the zeros, taking care to sum positive parts and negative parts first and then subtract in the end.
A second approach might be the double exponential quadrature formula adapted for oscillatory integrals. Please refer to Ooura and Mori, "The double exponential formula for oscillatory functions over the half infinite interval", J. Comp. App. Math. , 38, p353-360, 1991. Ooura has FORTRAN and C implementations on his website.
Finally, since you already know the correct answer to your integral, why do it numerically? Just to understand why
integral fails or is there a bigger picture?