The problem is to calculate the array of functions f(x) defined as:
Some coefficients of the matrix “A” are in the form:
For this reason, for large values “x” the matrix “A” becomes singular and the linear system bad conditioned.
Anyway, in double precision if “x > 706” the calculation of the coefficient gives underflow.
Matrix preconditioning seems not effective.
Are there math trick to overcome the problem?
Have you tried to obtain the LU decomposition and next use solve_triangular?
You can either do it manually with scipy.linalg.lu and scipy.linalg.solve_triangular or use the function dgesv. The first option will allow you to control better the middle steps, the second option is more confortable.
Note that if the matrix is defined positive, you can use the cholesky factorization instead of the LU.