During initial construction a dense matrix M is used. It seems that csr_matrix has been imported from scipy but is never actually used to instantiate a sparse csr matrix. Also, it seems that for sparse solving we are supplying the vector b as an array type. If this is the case, I believe scipy.sparse.linalg.spsolve will convert between dense to sparse and then back to dense during each time step because spsolve cannot assume that the solution X will be sparse.
I think to get the full power of scipy's solvers we might need to add a few changes here and during inital matrix M construction.
https://github.com/SimVascular/svZeroDSolver/blob/817b95f6aedde24f158b1e7148920e2cf912030f/svzerodsolver/time_integration.py#L63-L68
https://github.com/SimVascular/svZeroDSolver/blob/817b95f6aedde24f158b1e7148920e2cf912030f/svzerodsolver/time_integration.py#L93-L104
Avoiding excessive matrix reformating should significantly speed up sparse applications.
Let me know what y'all think @ktbolt @mrp089 @JonathanPham @richterjakob