A general program is being developed for the systematic and careful evaluation of Green functions in pion-nucleon problems. This paper is the first in a series of projects involving integral method for the calculation of the pion and nucleon propagators. A method is presented for evaluating efficiently and accurately integrals of the form integral r^m chi _l1(k₁r) chi _l2(k₂r). . . chi _ln(k_nr) dr, where m and n are arbitrary integers and chi _l(x) can be either a spherical Bessel function, j_l(x), or a spherical Neumann function, n_l(x). The range of integration depends upon the particular problem encountered. The prototype integral studied is I_ll'L(kk'P) identical to integral ₀^infinityr²j_l(kr) j_l'(k'r)j_L(pr)dr whose integrand for large r has a slowly decreasing oscillatory behaviour. Rapid convergence is ensured by rotating, in the complex plane, the upper part of this integral, giving an integrand which decreases exponentially. A scaling formula is used to evaluate I_ll'L(kk'p) for very small or very large values of the momenta. Also, it is shown that, if l, l' and L satisfy triangular inequalities and if l+l'+L is even, then k, k' and p must also satisfy triangular inequalities, which is the condition required by the vector delta function delta (k+k'-p). Finally the authors present sum rules and integral relations for I_ll'L(kk'p).