I want to solve linear hyperbolic system using Chebyshev collocation method. As this method puts severe constraint on the time step for the explicit time integration, I decided to switch to implicit time integration (backward Euler). I use time integrator
vode with options
order=1 from Python package
scipy.integrate. I do not supply Jacobian to the time integrator, only the vector of right-hand side. In the beginning it was great: the computation with resolution $N=100$ gave me solution 108 times faster with implicit time integration, then with explicit (Runge-Kutta 4th order).
Unfortunately, It turned out that if I increase the grid resolution, then computations become unstable—in initial time period solution is stable, but then it starts to "tremble" and develops instability.
This figure shows the numerically stable solution obtained with resolution $N=200$:
This figure shows the numerically unstable solution obtained with resolution $N=400$:
Could you please recommend to me implicit time integrators that are good to use with Chebyshev collocation method?