I have the following expression to be numerically integrated in a vector-based library (e.g. numpy, MATLAB, etc), $$ F(r_2) = \frac{1}{r_2^n} \int_0^{r_2} r_1^n f(r_1)\ \mathrm{d}r_1, $$ where $n$ is an integer that can have value up to $40$ and $r_1$ and $r_2$ can have values of many orders of magnitude (e.g. $r_2 \in [10^{-10},10^4]$).
The problem is, the numerical integration method that I use introduces a slight numerical instability which makes the integration result $\lim_{r_2\rightarrow 0}\int_0^{r_2}r_1^n f(r_1)\mathrm{d}r_1$ not exactly $0$, but something around $10^{-13}$. So when the integration result is divided by $r_2^n\sim10^{-10n}$ where $n \geq 3$, the result of $F(r_2\rightarrow 0)$ becomes really big.
I am looking for a way to allow the multiplication of $r_1^n/r_2^n$ takes place before the integration, so the numerical instability does not grow, or anything that prevents the numerical instability from the integration from growing really big.
The question is: how can I integrate the expression above numerically using numpy, MATLAB, or any vector-based library without growing the numerical instability?
P.S: I avoid for-loops if possible because I program in Python and for-loops can be really slow.