example.jl

using OrdinaryDiffEq, DiffEqOperators,Sundials

n = 600
M = rand(n,n)
K = rand(n,n)
ϵ = 1e-2
p = (ϵ,M,K)

u0 = ones(n)

function compute_residual!(resid,du,u,p,t)
        ϵ = p[1]
        M = p[2]
        K = p[3]
        resid .= M*du - ϵ*K*u
end

differential_vars = [true for i in 1:size(M)[1]]
prob = DAEProblem(compute_residual!,ϵ*M\(K*u0),u0,(0.0,1.0),differential_vars=differential_vars,p)
u = solve(prob,IDA(linear_solver=:GMRES))