Hey ya’ll, for the coupled channels scattering RBM stuff we’re having some difficulties using scipy.integrate for a coupled second order linear system, where each sub-matrix in the block matrix is sparse. (This corresponds to the Schroedinger eqn w/ coupling between channels).
I think the issue is that scipy is set up for first order systems, so for an ND 2nd order problem, you have to convert it to a 2ND 1st order problem. Beyond that, it just uses Runge-Kutta and iterates, which doesn’t take advantage of that fact that the system is linear and the sub matrices are sparse.
The end result is we’re getting some weird numerical instability and oscillations in the resulting wavefunctions, This could definitely be user error, but if anyone has experience with something similar I’d love to hear ideas!
I’ve been looking into using numba as well, so if anyone has experience with that too that would be great.