Assuming that you do not have to write the routine yourself, see the documentation for your preferred scientific subroutine package, write a calling program and test it out. There are any number of plotting and graphing packages that will display your data.
See “adaptive quadrature” in Wikipedia or “Numerical Recipes” by Press et al.
The basic idea is to apply a fixed N point quadrature rule (numerical integration) over a larger number of sub intervals until the results converge. Whole range, halves, quarters etc.
The last time I did this, about 10 years ago in Fortran, computational noise took over at 1024 intervals (or so).