The integrand is not singular at ##x = 0## so you could use a standard routine for adaptive quadrature, as implemented in the standard software (Maple, Mathematica, MATLAB, Octave, etc.). This will select a smaller mesh width near ##x = 0## and allow for a larger mesh width for large ##x##.
If you insist on doing it yourself, it is not so hard to implement one of the usual rules (e.g. trapezoidal or Simpson's) adaptively through recursion. Most NA books discuss this. If you need references, let me know. To deal with the limit at infinity, just restrict to a large but finite interval. Because of your behavior for large ##x##, it will be easy to get a bound on the contribution to the error.