The important functions in my program are "matrixHN" and "matrixSN". I use the numerical integration method "NewtonCotesRule", because it gives the best results. Is there a better integration method?
It would help if you provided some examples of the integrals you are using and what you may have already tried in Mathematica.