I have already tried using expoential solvers from ExponentialUtilities.jl , KrylovKit.jl and even solvers such as Exprb32 and 43 from DifferentialEquations.jl package for a miniature version of the original problem.
Here is a detailed reference to the same: Trying to solve a large system of highly stiff differential equations - #11 by balaji_sriram.
They are horribly slow in my case (and sometimes are highly inaccurate over a large time period- see the note in the next paragraph). Ive also tried various other stiff solvers previously (such as Rodas, Kencarp, Rosenbrock etc.), but found CVODE_BDF to be by far the most stable, accurate and fast solver of all. Since CVODE works well for my case, I would like to employ parallel processing (in my HPC) to these solvers (preferably). Also my HPC doesn’t have GPUs so at the moment I might not be able to explore the same.
A quick note of the system: The max values in the jacobian are of the order 1e12, the minimum nonzero values are 1e3 and i want to solve the problem till a span of 1e-3 seconds (a large time-span considering the fastest dynamics (1e12 Hz)are 9 orders faster than the slowest dynamics).
Also @ChrisRackauckas, as you recommended I will try using solvers such as LinearExponential(), and other State-Independent Solvers and let you know the results soon.