I’m trying to solve a DAE (mass matrix form) for a 2D spectral discretization. The number of DOF is going to be large, and the Jacobian, which I compute explicitly, is not especially sparse.

I’d like to use a preconditioned Krylov solver operating in a matrix-free mode, as applying the Jacobian can be done much faster than the equivalent matrix-vector multiplication. But I can’t figure out how to do that within DifferentialEquations.jl.

Going by the docs, the linsolve option to the native Julia stiff solvers apparently expects to pass matrices around. In SUNDIALS, which does support “ATimes” functions instead of matrices for its iterative solvers, I can’t figure out how to set them up.

Does anyone have a working example of using matrix-free methods in this context?