Rodas5 and ModelingToolkit much slower than Sundials' IDA on large systems

Thanks for the comments, I guess IDA is the way to go. Is there a way to define a DAEProblem in MTK to use with Sundials (I can’t really find any documentation about this)? When I naively try

	prob = DAEProblem(DAEFunction(sys), ICs, tspan, param_values)
	
	@time sol = solve(prob, IDA(), reltol = 1e-10, abstol = 1e-6);

I get the error ERROR: LoadError: MethodError: no method matching size(::Tuple{Float64, Float64}).

Full Trace
ERROR: LoadError: MethodError: no method matching size(::Tuple{Float64, Float64})
Closest candidates are:
  size(::Tuple, ::Integer) at tuple.jl:27
  size(::Union{LinearAlgebra.QR, LinearAlgebra.QRCompactWY, LinearAlgebra.QRPivoted}) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/qr.jl:524
  size(::Union{LinearAlgebra.QR, LinearAlgebra.QRCompactWY, LinearAlgebra.QRPivoted}, ::Integer) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/qr.jl:523
  ...
Stacktrace:
 [1] (DAEProblem{true, duType, tType, isinplace, P, F, K, D} where {duType, tType, isinplace, P, F, K, D})(f::DAEFunction{true, ModelingToolkit.var"#f#160"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##arg#576"), Symbol("##arg#577"), Symbol("##arg#578"), :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x77493f90, 0x4499b563, 0xe171eec5, 0xf4ea53b4, 0x48cc079b)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##out#579"), Symbol("##arg#576"), Symbol("##arg#577"), Symbol("##arg#578"), :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x3967f6b7, 0xd19b559f, 0xc7f890a0, 0xde6739f5, 0x3c2fea5e)}}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Nothing}, du0::Vector{Float64}, u0::Tuple{Float64, Float64}, tspan::Vector{Float64}, p::SciMLBase.NullParameters; differential_vars::Nothing, kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ SciMLBase ~/.julia/packages/SciMLBase/9EjAY/src/problems/dae_problems.jl:18
 [2] (DAEProblem{true, duType, tType, isinplace, P, F, K, D} where {duType, tType, isinplace, P, F, K, D})(f::DAEFunction{true, ModelingToolkit.var"#f#160"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##arg#576"), Symbol("##arg#577"), Symbol("##arg#578"), :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x77493f90, 0x4499b563, 0xe171eec5, 0xf4ea53b4, 0x48cc079b)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##out#579"), Symbol("##arg#576"), Symbol("##arg#577"), Symbol("##arg#578"), :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x3967f6b7, 0xd19b559f, 0xc7f890a0, 0xde6739f5, 0x3c2fea5e)}}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Nothing}, du0::Vector{Float64}, u0::Tuple{Float64, Float64}, tspan::Vector{Float64}, p::SciMLBase.NullParameters)
   @ SciMLBase ~/.julia/packages/SciMLBase/9EjAY/src/problems/dae_problems.jl:18
 [3] DAEProblem(f::DAEFunction{true, ModelingToolkit.var"#f#160"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##arg#576"), Symbol("##arg#577"), Symbol("##arg#578"), :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x77493f90, 0x4499b563, 0xe171eec5, 0xf4ea53b4, 0x48cc079b)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##out#579"), Symbol("##arg#576"), Symbol("##arg#577"), Symbol("##arg#578"), :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x3967f6b7, 0xd19b559f, 0xc7f890a0, 0xde6739f5, 0x3c2fea5e)}}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Nothing}, du0::Vector{Float64}, u0::Tuple{Float64, Float64}, tspan::Vector{Float64}, p::SciMLBase.NullParameters; kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ SciMLBase ~/.julia/packages/SciMLBase/9EjAY/src/problems/dae_problems.jl:37
 [4] DAEProblem(f::DAEFunction{true, ModelingToolkit.var"#f#160"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##arg#576"), Symbol("##arg#577"), Symbol("##arg#578"), :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x77493f90, 0x4499b563, 0xe171eec5, 0xf4ea53b4, 0x48cc079b)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##out#579"), Symbol("##arg#576"), Symbol("##arg#577"), Symbol("##arg#578"), :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x3967f6b7, 0xd19b559f, 0xc7f890a0, 0xde6739f5, 0x3c2fea5e)}}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Nothing}, du0::Vector{Float64}, u0::Tuple{Float64, Float64}, tspan::Vector{Float64}, p::SciMLBase.NullParameters) (repeats 2 times)
   @ SciMLBase ~/.julia/packages/SciMLBase/9EjAY/src/problems/dae_problems.jl:37
 [5] solve_system(sys::ODESystem; tspan::Tuple{Float64, Float64}, param_values::Vector{Float64})
   @ Main ~/Desktop/Silicon_Processes/Numerics/Example.jl:67
 [6] solve_system(sys::ODESystem)
   @ Main ~/Desktop/Silicon_Processes/Numerics/Example.jl:61
 [7] top-level scope
   @ ~/Desktop/Silicon_Processes/Numerics/Example.jl:78
 [8] include(fname::String)
   @ Base.MainInclude ./client.jl:444
 [9] top-level scope
   @ REPL[25]:1
in expression starting at /home/brady/Desktop/Silicon_Processes/Numerics/Example.jl:78

If not, would following something like
https://discourse.julialang.org/t/error-evaluating-sparse-analytical-jacobian-from-modelingtoolkit-in-dae-system/45764/2
be best for optimizing the Jacobian calculation and exploiting the sparsity?