D_sicm.jl

using ModelingToolkit
using ModelingToolkit: Differential, expand_derivatives, Expression, Operation, simplify_constants

function D(f)
    (args...) -> begin
        syms = map(_->Variable(gensym())(), args)
        ex = f(syms...)
        ds = map(s->Differential(s), syms)
        ops = [expand_derivatives(ds[i](ex)) for i in eachindex(syms)]
        rules = Dict( collect(zip(syms, args)) )
        foreach(op->_replace!(op, rules), ops)
        return map(simplify_constants, ops)
    end
end

function _replace!(op::Union{Vector{Expression}, Operation}, rules::Dict)
    if op isa Operation
        _replace!(op.args, rules)
    elseif eltype(op) <: Expression
        for i in eachindex(op)
            op[i] isa Operation && _replace!(op[i], rules)
            haskey(rules, op[i]) || continue
            op[i] = rules[op[i]]
        end
    else
        return nothing
    end
end