I want to solve numerically a system of equations defined as
\begin{align} \sum_{i,j} f_{ij}(r) &= 1 \\ 0 &= A_{ijkl}f_{kl}(r) + \frac{1}{r^6}B_{ijkl}f_{kl}(r) \\ &+ f_{kl}(r)\left[\int_{r_0}^{\infty}\frac{1}{r^4}f_{km}(r)\mathrm{d}r\frac{C_{kmon}}{\sqrt{f_{kk}(\infty)}}\right.\\ &+ \left. \int_{r_0}^{\infty}\frac{1}{r^4}f_{lm}(r)\mathrm{d}r\frac{C_{lmon}}{\sqrt{f_{ll}(\infty)}}\right] \end{align}
where the second equation holds for all $i,j,r$ and $A,B,C$ are known. My current numerical approach is to discretize over $r$, assume a starting set of $f_{ij}$ indep. of $r$, calculate the value between brackets and solve the resulting system for all $r$ to give a next guess for $f_{ij}$. This has a large radius of convergence but is too inaccurate.
What is the best numerical method/library for solving these types of equations, i.e. systems of equations involving integrals?
Numerical methods for nonlinear system of equations are problem dependent and require lots of turning. Your approach is correct in theory, but there are lots of convergence issues. Maybe this is the reason why there is no general method/library.
There are some iteration methods in
Scipy, which you can try and tune. It seems that to supply the Jacobian matrix is not essential in other direct or general methods, so maybe you can try them as well.