diff --git a/Project.toml b/Project.toml
index 110118a5..878777ba 100644
--- a/Project.toml
+++ b/Project.toml
@@ -29,6 +29,7 @@ Arpack = "0.5.1 - 0.5.3"
 DiffEqBase = "6.162"
 DiffEqCallbacks = "2, 3, 4"
 DiffEqNoiseProcess = "5.23.0"
+DifferentiationInterface = "0.6.39"
 FFTW = "1"
 FiniteDiff = "2.27.0"
 ForwardDiff = "0.10, 1"
@@ -52,6 +53,7 @@ julia = "1.10"
 
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
+DifferentiationInterface = "a0c0ee7d-e4b9-4e03-894e-1c5f64a51d63"
 FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41"
 JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -61,4 +63,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TestItemRunner = "f8b46487-2199-4994-9208-9a1283c18c0a"
 
 [targets]
-test = ["Aqua", "FiniteDiff", "JET", "LinearAlgebra", "SparseArrays", "Random", "Test", "TestItemRunner"]
+test = ["Aqua", "DifferentiationInterface", "FiniteDiff", "JET", "LinearAlgebra", "SparseArrays", "Random", "Test", "TestItemRunner"]
diff --git a/src/master.jl b/src/master.jl
index 8d17cbcb..9e8d523f 100644
--- a/src/master.jl
+++ b/src/master.jl
@@ -14,7 +14,7 @@ function master_h(tspan, rho0::Operator, H::AbstractOperator, J;
     _check_const.(J)
     _check_const.(Jdagger)
     check_master(rho0, H, J, Jdagger, rates)
-    tspan, rho0 = _promote_time_and_state(rho0, H, J, tspan)
+    tspan, rho0 = _promote_time_and_state(rho0, H, J, rates, tspan)
     tmp = copy(rho0)
     dmaster_(t, rho, drho) = dmaster_h!(drho, H, J, Jdagger, rates, rho, tmp)
     integrate_master(tspan, dmaster_, rho0, fout; kwargs...)
@@ -42,7 +42,7 @@ function master_nh(tspan, rho0::Operator, Hnh::AbstractOperator, J;
     _check_const.(J)
     _check_const.(Jdagger)
     check_master(rho0, Hnh, J, Jdagger, rates)
-    tspan, rho0 = _promote_time_and_state(rho0, Hnh, J, tspan)
+    tspan, rho0 = _promote_time_and_state(rho0, Hnh, J, rates, tspan)
     tmp = copy(rho0)
     dmaster_(t, rho, drho) = dmaster_nh!(drho, Hnh, Hnhdagger, J, Jdagger, rates, rho, tmp)
     integrate_master(tspan, dmaster_, rho0, fout; kwargs...)
@@ -88,7 +88,7 @@ function master(tspan, rho0::Operator, H::AbstractOperator, J;
     _check_const(H)
     _check_const.(J)
     _check_const.(Jdagger)
-    tspan, rho0 = _promote_time_and_state(rho0, H, J, tspan)
+    tspan, rho0 = _promote_time_and_state(rho0, H, J, rates, tspan)
     isreducible = check_master(rho0, H, J, Jdagger, rates)
     if !isreducible
         tmp = copy(rho0)
@@ -172,14 +172,32 @@ function master_nh_dynamic(tspan, rho0::Operator, f;
                 rates=nothing,
                 fout=nothing,
                 kwargs...)
+    tspan, rho0 = _promote_time_and_state(rho0, f, tspan)
     tmp = copy(rho0)
     dmaster_(t, rho, drho) = dmaster_nh_dynamic!(drho, f, rates, rho, tmp, t)
     integrate_master(tspan, dmaster_, rho0, fout; kwargs...)
 end
-
+function master_nh_dynamic(tspan, rho0::Operator, Hnh::AbstractTimeDependentOperator, J::Vector{<:AbstractTimeDependentOperator};
+    kwargs...)
+    promoted_tspan, rho0 = _promote_time_and_state(rho0, Hnh, J, tspan)
+    if promoted_tspan !== tspan # promote H
+        promoted_Hnh = TimeDependentSum(Hnh.coefficients, Hnh.static_op.operators; init_time=first(promoted_tspan))
+        promoted_J = [TimeDependentSum(Ji.coefficients, Ji.static_op.operators; init_time=first(promoted_tspan)) for Ji in J]
+        f = master_nh_dynamic_function(promoted_Hnh, promoted_J)
+    else
+        f = master_nh_dynamic_function(Hnh, J)
+    end
+    master_nh_dynamic(tspan, rho0, f; kwargs...)
+end
 function master_nh_dynamic(tspan, rho0::Operator, Hnh::AbstractTimeDependentOperator, J;
     kwargs...)
-    f = master_nh_dynamic_function(Hnh, J)
+    promoted_tspan, rho0 = _promote_time_and_state(rho0, Hnh, J, tspan)
+    if promoted_tspan !== tspan # promote H
+        promoted_Hnh = TimeDependentSum(Hnh.coefficients, Hnh.static_op.operators; init_time=first(promoted_tspan))
+        f = master_nh_dynamic_function(promoted_Hnh, J)
+    else
+        f = master_nh_dynamic_function(Hnh, J)
+    end
     master_nh_dynamic(tspan, rho0, f; kwargs...)
 end
 
@@ -218,6 +236,7 @@ function master_dynamic(tspan, rho0::Operator, f;
                 rates=nothing,
                 fout=nothing,
                 kwargs...)
+    tspan, rho0 = _promote_time_and_state(rho0, f, tspan)
     tmp = copy(rho0)
     dmaster_ = let f = f, tmp = tmp
         dmaster_(t, rho, drho) = dmaster_h_dynamic!(drho, f, rates, rho, tmp, t)
@@ -225,10 +244,27 @@ function master_dynamic(tspan, rho0::Operator, f;
     integrate_master(tspan, dmaster_, rho0, fout; kwargs...)
 end
 
+function master_dynamic(tspan, rho0::Operator, H::AbstractTimeDependentOperator, J::Vector{<:AbstractTimeDependentOperator};
+    kwargs...)
+    promoted_tspan, rho0 = _promote_time_and_state(rho0, H, J, tspan)
+    if promoted_tspan !== tspan # promote H
+        promoted_H = TimeDependentSum(H.coefficients, H.static_op.operators; init_time=first(promoted_tspan))
+        promoted_J = [TimeDependentSum(Ji.coefficients, Ji.static_op.operators; init_time=first(promoted_tspan)) for Ji in J]
+        f = master_h_dynamic_function(promoted_H, promoted_J)
+    else
+        f = master_h_dynamic_function(H, J)
+    end
+    return master_dynamic(promoted_tspan, rho0, f; kwargs...)
+end
 function master_dynamic(tspan, rho0::Operator, H::AbstractTimeDependentOperator, J;
     kwargs...)
-    f = master_h_dynamic_function(H, J)
-    master_dynamic(tspan, rho0, f; kwargs...)
+    promoted_tspan, rho0 = _promote_time_and_state(rho0, H, J, tspan)
+    if promoted_tspan !== tspan # promote H
+        promoted_H = TimeDependentSum(H.coefficients, H.static_op.operators; init_time=first(promoted_tspan))
+        return master_dynamic(promoted_tspan, rho0, master_h_dynamic_function(promoted_H, J); kwargs...)
+    else
+        return master_dynamic(promoted_tspan, rho0, master_h_dynamic_function(H, J); kwargs...)
+    end
 end
 
 # Automatically convert Ket states to density operators
diff --git a/src/timeevolution_base.jl b/src/timeevolution_base.jl
index e90ecf99..dc2e8c2e 100644
--- a/src/timeevolution_base.jl
+++ b/src/timeevolution_base.jl
@@ -143,8 +143,29 @@ function _promote_time_and_state(u0, H::AbstractOperator, J, tspan)
     tspan_promote = DiffEqBase.promote_tspan(u0_promote.data, p, tspan, nothing, Dict{Symbol, Any}())
     return tspan_promote, u0_promote
 end
-
-_promote_time_and_state(u0, f, tspan) = _promote_time_and_state(u0, f(first(tspan)..., u0), tspan)
+function _promote_time_and_state(u0, H::AbstractOperator, J, rates, tspan)
+    # TODO: Find an alternative to promote_dual, which was moved to
+    #       an extension in DiffEqBase 6.162.0
+    ext = Base.get_extension(DiffEqBase, :DiffEqBaseForwardDiffExt)
+    Ts = reduce(ext.promote_dual, (eltype(H), DiffEqBase.anyeltypedual(J), DiffEqBase.anyeltypedual(rates)))
+    Tt = real(Ts)
+    p = Vector{Tt}(undef,0)
+    u0_promote = DiffEqBase.promote_u0(u0, p, tspan[1])
+    tspan_promote = DiffEqBase.promote_tspan(u0_promote.data, p, tspan, nothing, Dict{Symbol, Any}())
+    return tspan_promote, u0_promote
+end
+_promote_time_and_state(u0, f::Function, tspan) = _promote_time_and_state(u0, f(first(tspan)..., u0), tspan)
+function _promote_time_and_state(u0, f::Union{Tuple, Vector}, tspan)
+    # TODO: Find an alternative to promote_dual, which was moved to
+    #       an extension in DiffEqBase 6.162.0
+    ext = Base.get_extension(DiffEqBase, :DiffEqBaseForwardDiffExt)
+    Ts = reduce(ext.promote_dual, (eltype(f[1]), DiffEqBase.anyeltypedual.(f[2:end])...))
+    Tt = real(Ts)
+    p = Vector{Tt}(undef,0)
+    u0_promote = DiffEqBase.promote_u0(u0, p, tspan[1])
+    tspan_promote = DiffEqBase.promote_tspan(u0_promote.data, p, tspan, nothing, Dict{Symbol, Any}())
+    return tspan_promote, u0_promote
+end
 
 @inline function DiffEqBase.promote_u0(u0::Ket, p, t0)
     u0data_promote = DiffEqBase.promote_u0(u0.data, p, t0)
diff --git a/test/test_autodiff.jl b/test/test_autodiff.jl
new file mode 100644
index 00000000..ab3e0fd0
--- /dev/null
+++ b/test/test_autodiff.jl
@@ -0,0 +1,501 @@
+@testitem "Test automatic differentation backends with QuantumOptics.jl" begin
+    using Test
+    using QuantumOptics
+    using DifferentiationInterface
+    import ForwardDiff
+    import FiniteDiff
+    
+    # cost parameters
+    pnum, pvec = rand(Float64), rand(Float64, 2)
+    # in-place buffers
+    bufvec, bufmat = similar(pvec), similar(pvec * pvec')
+    # random system definitions
+    basis = SpinBasis(10//1)
+    ψ0 = randstate(basis)
+    ρ0 = dm(ψ0)
+    H1, H2 = randoperator(basis), randoperator(basis)
+    H1t, H2t = TimeDependentSum(sin=>H1), TimeDependentSum(sin=>H2)
+    tspan = [0.0, 1.0]
+
+    @testset "time-evolution" begin
+
+        @testset "schroedinger" begin
+
+            function cost(p::Real)
+                _, ψt = timeevolution.schroedinger(tspan, ψ0, p*H1)
+                return 1 - norm(ψt)^2
+            end
+            function costvec(p::Real)
+                _, ψt = timeevolution.schroedinger(tspan, ψ0, p*H1)
+                return [1 - norm(ψt)^2, norm(ψt)]
+            end
+            function cost(p::Vector)
+                _, ψt = timeevolution.schroedinger(tspan, ψ0, p[1]*H1 + p[2]*H2)
+                return 1 - norm(ψt)^2
+            end
+            function costvec(p::Vector)
+                _, ψt = timeevolution.schroedinger(tspan, ψ0, p[1]*H1 + p[2]*H2)
+                return [1 - norm(ψt)^2, norm(ψt)]
+            end
+            for backend in [AutoForwardDiff(), AutoFiniteDiff()]
+                # out-of-place
+                @test_nowarn derivative(cost, backend, pnum)
+                @test_nowarn second_derivative(cost, backend, pnum)
+                @test_nowarn gradient(cost, backend, pvec)
+                @test_nowarn jacobian(costvec, backend, pvec)
+                @test_nowarn hessian(cost, backend, pvec)
+                # in-place
+                @test_nowarn derivative!(costvec, bufvec, backend, pnum)
+                @test_nowarn second_derivative!(costvec, bufvec, backend, pnum)
+                @test_nowarn gradient!(cost, bufvec, backend, pvec)
+                @test_nowarn jacobian!(costvec, bufmat, backend, pvec)
+                @test_nowarn hessian!(cost, bufmat, backend, pvec)
+            end
+        end
+
+        @testset "dynamic schroedinger" begin
+
+            @testset "time-dependent functions" begin
+
+                function cost(p::Real)
+                    _, ψt = timeevolution.schroedinger_dynamic(tspan, ψ0, (t, ψ) -> exp(p*t)*H1)
+                    return 1 - norm(ψt)^2
+                end
+                function costvec(p::Real)
+                    _, ψt = timeevolution.schroedinger_dynamic(tspan, ψ0, (t, ψ) -> exp(p*t)*H1)
+                    return [1 - norm(ψt)^2, norm(ψt)]
+                end
+                function cost(p::Vector)
+                    _, ψt = timeevolution.schroedinger_dynamic(tspan, ψ0, (t, ψ) -> exp(p[1]*t)*H1 + exp(p[2]*t)*H2)
+                    return 1 - norm(ψt)^2
+                end
+                function costvec(p::Vector)
+                    _, ψt = timeevolution.schroedinger_dynamic(tspan, ψ0, (t, ψ) -> exp(p[1]*t)*H1 + exp(p[2]*t)*H2)
+                    return [1 - norm(ψt)^2, norm(ψt)]
+                end
+                for backend in [AutoForwardDiff(), AutoFiniteDiff()]
+                    # out-of-place
+                    @test_nowarn derivative(cost, backend, pnum)
+                    @test_nowarn second_derivative(cost, backend, pnum)
+                    @test_nowarn gradient(cost, backend, pvec)
+                    @test_nowarn jacobian(costvec, backend, pvec)
+                    @test_nowarn hessian(cost, backend, pvec)
+                    # in-place
+                    @test_nowarn derivative!(costvec, bufvec, backend, pnum)
+                    @test_nowarn second_derivative!(costvec, bufvec, backend, pnum)
+                    @test_nowarn gradient!(cost, bufvec, backend, pvec)
+                    @test_nowarn jacobian!(costvec, bufmat, backend, pvec)
+                    @test_nowarn hessian!(cost, bufmat, backend, pvec)
+                end
+            end
+
+            @testset "TimeDependentOperators" begin
+
+                function cost(p::Real)
+                    _, ψt = timeevolution.schroedinger_dynamic(tspan, ψ0, p*H1t)
+                    return 1 - norm(ψt)^2
+                end
+                function costvec(p::Real)
+                    _, ψt = timeevolution.schroedinger_dynamic(tspan, ψ0, p*H1t)
+                    return [1 - norm(ψt)^2, norm(ψt)]
+                end
+                function cost(p::Vector)
+                    _, ψt = timeevolution.schroedinger_dynamic(tspan, ψ0, p[1]*H1t + p[2]*H2t)
+                    return 1 - norm(ψt)^2
+                end
+                function costvec(p::Vector)
+                    _, ψt = timeevolution.schroedinger_dynamic(tspan, ψ0, p[1]*H1t + p[2]*H2t)
+                    return [1 - norm(ψt)^2, norm(ψt)]
+                end
+                for backend in [AutoForwardDiff(), AutoFiniteDiff()]
+                    # out-of-place
+                    @test_nowarn derivative(cost, backend, pnum)
+                    @test_nowarn second_derivative(cost, backend, pnum)
+                    @test_nowarn gradient(cost, backend, pvec)
+                    @test_nowarn jacobian(costvec, backend, pvec)
+                    @test_nowarn hessian(cost, backend, pvec)
+                    # in-place
+                    @test_nowarn derivative!(costvec, bufvec, backend, pnum)
+                    @test_nowarn second_derivative!(costvec, bufvec, backend, pnum)
+                    @test_nowarn gradient!(cost, bufvec, backend, pvec)
+                    @test_nowarn jacobian!(costvec, bufmat, backend, pvec)
+                    @test_nowarn hessian!(cost, bufmat, backend, pvec)
+                end
+            end
+        end
+
+        for eq in [:master, :master_h, :master_nh]
+
+            @testset "$eq" begin
+
+                # differentiate Hamiltonian
+                @testset "Hamiltonian" begin
+
+                    function cost(p::Real)
+                         _, ρt = @eval(timeevolution.$eq)(tspan, ρ0, p*H1, [])
+                        return 1 - norm(ρt)^2
+                    end
+                    function costvec(p::Real)
+                        _, ψt = @eval(timeevolution.$eq)(tspan, ρ0, p*H1, [])
+                        return [1 - norm(ψt)^2, norm(ψt)]
+                    end
+                    function cost(p::Vector)
+                        _, ρt = @eval(timeevolution.$eq)(tspan, ρ0, p[1]*H1 + p[2]*H2, [])
+                        return 1 - norm(ρt)^2
+                    end
+                    function costvec(p::Vector)
+                        _, ρt = @eval(timeevolution.$eq)(tspan, ρ0, p[1]*H1 + p[2]*H2, [])
+                        return [1 - norm(ρt)^2, norm(ρt)]
+                    end
+                    for backend in [AutoForwardDiff(), AutoFiniteDiff()]
+                        # out-of-place
+                        @test_nowarn derivative(cost, backend, pnum)
+                        @test_nowarn second_derivative(cost, backend, pnum)
+                        @test_nowarn gradient(cost, backend, pvec)
+                        @test_nowarn jacobian(costvec, backend, pvec)
+                        @test_nowarn hessian(cost, backend, pvec)
+                        # in-place
+                        @test_nowarn derivative!(costvec, bufvec, backend, pnum)
+                        @test_nowarn second_derivative!(costvec, bufvec, backend, pnum)
+                        @test_nowarn gradient!(cost, bufvec, backend, pvec)
+                        @test_nowarn jacobian!(costvec, bufmat, backend, pvec)
+                        @test_nowarn hessian!(cost, bufmat, backend, pvec)
+                    end
+                end
+
+                # differentiate jump operators
+                @testset "jump operators" begin
+
+                    function cost(p::Real)
+                        _, ρt = @eval(timeevolution.$eq)(tspan, ρ0, H1, [p*H2])
+                        return 1 - norm(ρt)^2
+                    end
+                    function costvec(p::Real)
+                        _, ψt = @eval(timeevolution.$eq)(tspan, ρ0, H1, [p*H2])
+                        return [1 - norm(ψt)^2, norm(ψt)]
+                    end
+                    function cost(p::Vector)
+                        _, ρt = @eval(timeevolution.$eq)(tspan, ρ0, H1, [p[1]*H1, p[2]*H2])
+                        return 1 - norm(ρt)^2
+                    end
+                    function costvec(p::Vector)
+                        _, ρt = @eval(timeevolution.$eq)(tspan, ρ0, H1, [p[1]*H1, p[2]*H2])
+                        return [1 - norm(ρt)^2, norm(ρt)]
+                    end
+                    for backend in [AutoForwardDiff(), AutoFiniteDiff()]
+                        # out-of-place
+                        @test_nowarn derivative(cost, backend, pnum)
+                        @test_nowarn second_derivative(cost, backend, pnum)
+                        @test_nowarn gradient(cost, backend, pvec)
+                        @test_nowarn jacobian(costvec, backend, pvec)
+                        @test_nowarn hessian(cost, backend, pvec)
+                        # in-place
+                        @test_nowarn derivative!(costvec, bufvec, backend, pnum)
+                        @test_nowarn second_derivative!(costvec, bufvec, backend, pnum)
+                        @test_nowarn gradient!(cost, bufvec, backend, pvec)
+                        @test_nowarn jacobian!(costvec, bufmat, backend, pvec)
+                        @test_nowarn hessian!(cost, bufmat, backend, pvec)
+                    end
+                end
+            end
+        end
+
+        @testset "master_dynamic" begin
+
+            @testset "time-dependent functions" begin
+
+                # differentiate Hamiltonian
+                @testset "Hamiltonian" begin
+
+                    function cost(p::Real)
+                        _, ρt = timeevolution.master_dynamic(tspan, ρ0, (t, ρ) -> (exp(p*t)*H1, [], []))
+                        return 1 - norm(ρt)^2
+                    end
+                    function costvec(p::Real)
+                        _, ψt = timeevolution.master_dynamic(tspan, ρ0, (t, ρ) -> (exp(p*t)*H1, [], []))
+                        return [1 - norm(ψt)^2, norm(ψt)]
+                    end
+                    function cost(p::Vector)
+                        _, ρt = timeevolution.master_dynamic(tspan, ρ0, (t, ρ) -> (exp(p[1]*t)*H1 + exp(p[2]*t)*H2, [], []))
+                        return 1 - norm(ρt)^2
+                    end
+                    function costvec(p::Vector)
+                        _, ρt = timeevolution.master_dynamic(tspan, ρ0, (t, ρ) -> (exp(p[1]*t)*H1 + exp(p[2]*t)*H2, [], []))
+                        return [1 - norm(ρt)^2, norm(ρt)]
+                    end
+                    for backend in [AutoForwardDiff(), AutoFiniteDiff()]
+                        # out-of-place
+                        @test_nowarn derivative(cost, backend, pnum)
+                        @test_nowarn second_derivative(cost, backend, pnum)
+                        @test_nowarn gradient(cost, backend, pvec)
+                        @test_nowarn jacobian(costvec, backend, pvec)
+                        @test_nowarn hessian(cost, backend, pvec)
+                        # in-place
+                        @test_nowarn derivative!(costvec, bufvec, backend, pnum)
+                        @test_nowarn second_derivative!(costvec, bufvec, backend, pnum)
+                        @test_nowarn gradient!(cost, bufvec, backend, pvec)
+                        @test_nowarn jacobian!(costvec, bufmat, backend, pvec)
+                        @test_nowarn hessian!(cost, bufmat, backend, pvec)
+                    end
+                end
+
+                # differentiate jump operators
+                @testset "jump operators" begin
+
+                    function cost(p::Real)
+                        _, ρt = timeevolution.master_dynamic(tspan, ρ0, (t, ρ) -> (H1, [exp(p*t)*H2], [dagger(exp(p*t)*H2)]))
+                        return 1 - norm(ρt)^2
+                    end
+                    function costvec(p::Real)
+                        _, ψt = timeevolution.master_dynamic(tspan, ρ0, (t, ρ) -> (H1, [exp(p*t)*H2], [dagger(exp(p*t)*H2)]))
+                        return [1 - norm(ψt)^2, norm(ψt)]
+                    end
+                    function cost(p::Vector)
+                        _, ρt = timeevolution.master_dynamic(tspan, ρ0, (t, ρ) -> (H1, [exp(p[1]*t)*H1, exp(p[2]*t)*H2], [dagger(exp(p[1]*t)*H1), dagger(exp(p[2]*t)*H2)]))
+                        return 1 - norm(ρt)^2
+                    end
+                    function costvec(p::Vector)
+                        _, ρt = timeevolution.master_dynamic(tspan, ρ0, (t, ρ) -> (H1, [exp(p[1]*t)*H1, exp(p[2]*t)*H2], [dagger(exp(p[1]*t)*H1), dagger(exp(p[2]*t)*H2)]))
+                        return [1 - norm(ρt)^2, norm(ρt)]
+                    end
+                    for backend in [AutoForwardDiff(), AutoFiniteDiff()]
+                        # out-of-place
+                        @test_nowarn derivative(cost, backend, pnum)
+                        @test_nowarn second_derivative(cost, backend, pnum)
+                        @test_nowarn gradient(cost, backend, pvec)
+                        @test_nowarn jacobian(costvec, backend, pvec)
+                        @test_nowarn hessian(cost, backend, pvec)
+                        # in-place
+                        @test_nowarn derivative!(costvec, bufvec, backend, pnum)
+                        @test_nowarn second_derivative!(costvec, bufvec, backend, pnum)
+                        @test_nowarn gradient!(cost, bufvec, backend, pvec)
+                        @test_nowarn jacobian!(costvec, bufmat, backend, pvec)
+                        @test_nowarn hessian!(cost, bufmat, backend, pvec)
+                    end
+                end
+
+                # differentiate rates
+                @testset "rates" begin
+
+                    function cost(p::Real)
+                        _, ρt = timeevolution.master_dynamic(tspan, ρ0, (t, ρ) -> (H1, [H2], [dagger(H2)], [exp(p*t)]))
+                        return 1 - norm(ρt)^2
+                    end
+                    function costvec(p::Real)
+                        _, ψt = timeevolution.master_dynamic(tspan, ρ0, (t, ρ) -> (H1, [H2], [dagger(H2)], [exp(p*t)]))
+                        return [1 - norm(ψt)^2, norm(ψt)]
+                    end
+                    function cost(p::Vector)
+                        _, ρt = timeevolution.master_dynamic(tspan, ρ0, (t, ρ) -> (H1, [H1, H2], [dagger(H1), dagger(H2)], [exp(p[1]*t), exp(p[2]*t)]))
+                        return 1 - norm(ρt)^2
+                    end
+                    function costvec(p::Vector)
+                        _, ρt = timeevolution.master_dynamic(tspan, ρ0, (t, ρ) -> (H1, [H1, H2], [dagger(H1), dagger(H2)], [exp(p[1]*t), exp(p[2]*t)]))
+                        return [1 - norm(ρt)^2, norm(ρt)]
+                    end
+                    for backend in [AutoForwardDiff(), AutoFiniteDiff()]
+                        # out-of-place
+                        @test_nowarn derivative(cost, backend, pnum)
+                        @test_nowarn second_derivative(cost, backend, pnum)
+                        @test_nowarn gradient(cost, backend, pvec)
+                        @test_nowarn jacobian(costvec, backend, pvec)
+                        @test_nowarn hessian(cost, backend, pvec)
+                        # in-place
+                        @test_nowarn derivative!(costvec, bufvec, backend, pnum)
+                        @test_nowarn second_derivative!(costvec, bufvec, backend, pnum)
+                        @test_nowarn gradient!(cost, bufvec, backend, pvec)
+                        @test_nowarn jacobian!(costvec, bufmat, backend, pvec)
+                        @test_nowarn hessian!(cost, bufmat, backend, pvec)
+                    end
+                end
+            end
+
+            @testset "TimeDependentOperators" begin
+
+                # differentiate Hamiltonian
+                @testset "Hamiltonian" begin
+
+                    function cost(p::Real)
+                        _, ρt = timeevolution.master_dynamic(tspan, ρ0, p*H1t, [p*H1t])
+                        return 1 - norm(ρt)^2
+                    end
+                    function costvec(p::Real)
+                        _, ψt = timeevolution.master_dynamic(tspan, ρ0, p*H1t, [p*H1t])
+                        return [1 - norm(ψt)^2, norm(ψt)]
+                    end
+                    function cost(p::Vector)
+                        _, ρt = timeevolution.master_dynamic(tspan, ρ0, p[1]*H1t, [p[2]*H1t])
+                        return 1 - norm(ρt)^2
+                    end
+                    function costvec(p::Vector)
+                        _, ρt = timeevolution.master_dynamic(tspan, ρ0, p[1]*H1t, [p[1]*H1t])
+                        return [1 - norm(ρt)^2, norm(ρt)]
+                    end
+                    for backend in [AutoForwardDiff(), AutoFiniteDiff()]
+                        # out-of-place
+                        @test_nowarn derivative(cost, backend, pnum)
+                        @test_nowarn second_derivative(cost, backend, pnum)
+                        @test_nowarn gradient(cost, backend, pvec)
+                        @test_nowarn jacobian(costvec, backend, pvec)
+                        @test_nowarn hessian(cost, backend, pvec)
+                        # in-place
+                        @test_nowarn derivative!(costvec, bufvec, backend, pnum)
+                        @test_nowarn second_derivative!(costvec, bufvec, backend, pnum)
+                        @test_nowarn gradient!(cost, bufvec, backend, pvec)
+                        @test_nowarn jacobian!(costvec, bufmat, backend, pvec)
+                        @test_nowarn hessian!(cost, bufmat, backend, pvec)
+                    end
+                end
+            end
+        end
+
+        @testset "master_nh_dynamic" begin
+
+            @testset "time-dependent functions" begin
+
+                # differentiate Hamiltonian
+                @testset "Hamiltonian" begin
+
+                    function cost(p::Real)
+                        _, ρt = timeevolution.master_nh_dynamic(tspan, ρ0, (t, ρ) -> (exp(p*t)*H1, dagger(exp(p*t)*H1), [], []))
+                        return 1 - norm(ρt)^2
+                    end
+                    function costvec(p::Real)
+                        _, ψt = timeevolution.master_nh_dynamic(tspan, ρ0, (t, ρ) -> (exp(p*t)*H1, dagger(exp(p*t)*H1), [], []))
+                        return [1 - norm(ψt)^2, norm(ψt)]
+                    end
+                    function cost(p::Vector)
+                        _, ρt = timeevolution.master_nh_dynamic(tspan, ρ0, (t, ρ) -> (exp(p[1]*t)*H1 + exp(p[2]*t)*H2, dagger(exp(p[1]*t)*H1 + exp(p[2]*t)*H2), [], []))
+                        return 1 - norm(ρt)^2
+                    end
+                    function costvec(p::Vector)
+                        _, ρt = timeevolution.master_nh_dynamic(tspan, ρ0, (t, ρ) -> (exp(p[1]*t)*H1 + exp(p[2]*t)*H2, dagger(exp(p[1]*t)*H1 + exp(p[2]*t)*H2), [], []))
+                        return [1 - norm(ρt)^2, norm(ρt)]
+                    end
+                    for backend in [AutoForwardDiff(), AutoFiniteDiff()]
+                        # out-of-place
+                        @test_nowarn derivative(cost, backend, pnum)
+                        @test_nowarn second_derivative(cost, backend, pnum)
+                        @test_nowarn gradient(cost, backend, pvec)
+                        @test_nowarn jacobian(costvec, backend, pvec)
+                        @test_nowarn hessian(cost, backend, pvec)
+                        # in-place
+                        @test_nowarn derivative!(costvec, bufvec, backend, pnum)
+                        @test_nowarn second_derivative!(costvec, bufvec, backend, pnum)
+                        @test_nowarn gradient!(cost, bufvec, backend, pvec)
+                        @test_nowarn jacobian!(costvec, bufmat, backend, pvec)
+                        @test_nowarn hessian!(cost, bufmat, backend, pvec)
+                    end
+                end
+
+                # differentiate jump operators
+                @testset "jump operators" begin
+
+                    function cost(p::Real)
+                        _, ρt = timeevolution.master_nh_dynamic(tspan, ρ0, (t, ρ) -> (H1, dagger(H1), [exp(p*t)*H2], [dagger(exp(p*t)*H2)]))
+                        return 1 - norm(ρt)^2
+                    end
+                    function costvec(p::Real)
+                        _, ψt = timeevolution.master_nh_dynamic(tspan, ρ0, (t, ρ) -> (H1, dagger(H1), [exp(p*t)*H2], [dagger(exp(p*t)*H2)]))
+                        return [1 - norm(ψt)^2, norm(ψt)]
+                    end
+                    function cost(p::Vector)
+                        _, ρt = timeevolution.master_nh_dynamic(tspan, ρ0, (t, ρ) -> (H1, dagger(H1), [exp(p[1]*t)*H1, exp(p[2]*t)*H2], [dagger(exp(p[1]*t)*H1), dagger(exp(p[2]*t)*H2)]))
+                        return 1 - norm(ρt)^2
+                    end
+                    function costvec(p::Vector)
+                        _, ρt = timeevolution.master_nh_dynamic(tspan, ρ0, (t, ρ) -> (H1, dagger(H1), [exp(p[1]*t)*H1, exp(p[2]*t)*H2], [dagger(exp(p[1]*t)*H1), dagger(exp(p[2]*t)*H2)]))
+                        return [1 - norm(ρt)^2, norm(ρt)]
+                    end
+                    for backend in [AutoForwardDiff(), AutoFiniteDiff()]
+                        # out-of-place
+                        @test_nowarn derivative(cost, backend, pnum)
+                        @test_nowarn second_derivative(cost, backend, pnum)
+                        @test_nowarn gradient(cost, backend, pvec)
+                        @test_nowarn jacobian(costvec, backend, pvec)
+                        @test_nowarn hessian(cost, backend, pvec)
+                        # in-place
+                        @test_nowarn derivative!(costvec, bufvec, backend, pnum)
+                        @test_nowarn second_derivative!(costvec, bufvec, backend, pnum)
+                        @test_nowarn gradient!(cost, bufvec, backend, pvec)
+                        @test_nowarn jacobian!(costvec, bufmat, backend, pvec)
+                        @test_nowarn hessian!(cost, bufmat, backend, pvec)
+                    end
+                end
+
+                # differentiate rates
+                @testset "rates" begin
+
+                    function cost(p::Real)
+                        _, ρt = timeevolution.master_nh_dynamic(tspan, ρ0, (t, ρ) -> (H1, dagger(H1), [H2], [dagger(H2)], [exp(p*t)]))
+                        return 1 - norm(ρt)^2
+                    end
+                    function costvec(p::Real)
+                        _, ψt = timeevolution.master_nh_dynamic(tspan, ρ0, (t, ρ) -> (H1, dagger(H1), [H2], [dagger(H2)], [exp(p*t)]))
+                        return [1 - norm(ψt)^2, norm(ψt)]
+                    end
+                    function cost(p::Vector)
+                        _, ρt = timeevolution.master_nh_dynamic(tspan, ρ0, (t, ρ) -> (H1, dagger(H1), [H1, H2], [dagger(H1), dagger(H2)], [exp(p[1]*t), exp(p[2]*t)]))
+                        return 1 - norm(ρt)^2
+                    end
+                    function costvec(p::Vector)
+                        _, ρt = timeevolution.master_nh_dynamic(tspan, ρ0, (t, ρ) -> (H1, dagger(H1), [H1, H2], [dagger(H1), dagger(H2)], [exp(p[1]*t), exp(p[2]*t)]))
+                        return [1 - norm(ρt)^2, norm(ρt)]
+                    end
+                    for backend in [AutoForwardDiff(), AutoFiniteDiff()]
+                        # out-of-place
+                        @test_nowarn derivative(cost, backend, pnum)
+                        @test_nowarn second_derivative(cost, backend, pnum)
+                        @test_nowarn gradient(cost, backend, pvec)
+                        @test_nowarn jacobian(costvec, backend, pvec)
+                        @test_nowarn hessian(cost, backend, pvec)
+                        # in-place
+                        @test_nowarn derivative!(costvec, bufvec, backend, pnum)
+                        @test_nowarn second_derivative!(costvec, bufvec, backend, pnum)
+                        @test_nowarn gradient!(cost, bufvec, backend, pvec)
+                        @test_nowarn jacobian!(costvec, bufmat, backend, pvec)
+                        @test_nowarn hessian!(cost, bufmat, backend, pvec)
+                    end
+                end
+            end
+
+            @testset "TimeDependentOperators" begin
+
+                # differentiate Hamiltonian
+                @testset "Hamiltonian" begin
+
+                    function cost(p::Real)
+                        _, ρt = timeevolution.master_nh_dynamic(tspan, ρ0, p*H1t, [p*H1t])
+                        return 1 - norm(ρt)^2
+                    end
+                    function costvec(p::Real)
+                        _, ψt = timeevolution.master_nh_dynamic(tspan, ρ0, p*H1t, [p*H1t])
+                        return [1 - norm(ψt)^2, norm(ψt)]
+                    end
+                    function cost(p::Vector)
+                        _, ρt = timeevolution.master_nh_dynamic(tspan, ρ0, p[1]*H1t, [p[2]*H1t])
+                        return 1 - norm(ρt)^2
+                    end
+                    function costvec(p::Vector)
+                        _, ρt = timeevolution.master_nh_dynamic(tspan, ρ0, p[1]*H1t, [p[1]*H1t])
+                        return [1 - norm(ρt)^2, norm(ρt)]
+                    end
+                    for backend in [AutoForwardDiff(), AutoFiniteDiff()]
+                        # out-of-place
+                        @test_nowarn derivative(cost, backend, pnum)
+                        @test_nowarn second_derivative(cost, backend, pnum)
+                        @test_nowarn gradient(cost, backend, pvec)
+                        @test_nowarn jacobian(costvec, backend, pvec)
+                        @test_nowarn hessian(cost, backend, pvec)
+                        # in-place
+                        @test_nowarn derivative!(costvec, bufvec, backend, pnum)
+                        @test_nowarn second_derivative!(costvec, bufvec, backend, pnum)
+                        @test_nowarn gradient!(cost, bufvec, backend, pvec)
+                        @test_nowarn jacobian!(costvec, bufmat, backend, pvec)
+                        @test_nowarn hessian!(cost, bufmat, backend, pvec)
+                    end
+                end
+            end
+        end
+    end
+end
\ No newline at end of file