fn solve_interacting( local: LocalMatrix, coupling: DMatrix<Complex>, ) -> Result<InteractionResidual>