@testitem "Double Pendulum" begin
    using ForwardDiff, Random

    Random.seed!(1)
    T = Float64
    system = StabilizedNDEs.DoublePendulum{T}(rand(T, 4)...)

    # One test for each experiment version
    u = rand(T, 8)
    jac_forward_diff_1 =
        ForwardDiff.jacobian(u -> StabilizedNDEs.invariants(u, nothing, system, Val(1)), u)
    jac_analytic_1 = StabilizedNDEs.invariants_jacobian(u, nothing, system, Val(1))
    @test jac_analytic_1 ≈ jac_forward_diff_1 atol = 1e-6

    u = rand(T, 8)
    jac_forward_diff_2 =
        ForwardDiff.jacobian(u -> StabilizedNDEs.invariants(u, nothing, system, Val(2)), u)
    jac_analytic_2 = StabilizedNDEs.invariants_jacobian(u, nothing, system, Val(2))
    @test jac_analytic_2 ≈ jac_forward_diff_2 atol = 1e-6

    u = rand(T, 4)
    jac_forward_diff_3 =
        ForwardDiff.jacobian(u -> StabilizedNDEs.invariants(u, nothing, system, Val(3)), u)
    jac_analytic_3 = StabilizedNDEs.invariants_jacobian(u, nothing, system, Val(3))
    @test jac_analytic_3 ≈ jac_forward_diff_3 atol = 1e-6

    u = rand(T, 8)
    jac_forward_diff_4 =
        ForwardDiff.jacobian(u -> StabilizedNDEs.invariants(u, nothing, system, Val(4)), u)
    jac_analytic_4 = StabilizedNDEs.invariants_jacobian(u, nothing, system, Val(4))
    @test jac_analytic_4 ≈ jac_forward_diff_4 atol = 1e-6

    u = rand(T, 8)
    jac_forward_diff_5 =
        ForwardDiff.jacobian(u -> StabilizedNDEs.invariants(u, nothing, system, Val(5)), u)
    jac_analytic_5 = StabilizedNDEs.invariants_jacobian(u, nothing, system, Val(5))
    @test jac_analytic_5 ≈ jac_forward_diff_5 atol = 1e-6
end

@testitem "Kepler Problem" begin
    using ForwardDiff, Random

    Random.seed!(1)
    T = Float64
    system = StabilizedNDEs.KeplerProblem{T}()

    # One test for each experiment version
    u = rand(T, 4)
    jac_forward_diff_1 =
        ForwardDiff.jacobian(u -> StabilizedNDEs.invariants(u, nothing, system, Val(1)), u)
    jac_analytic_1 = StabilizedNDEs.invariants_jacobian(u, nothing, system, Val(1))
    @test jac_analytic_1 ≈ jac_forward_diff_1 atol = 1e-6

    u = rand(T, 4)
    jac_forward_diff_2 =
        ForwardDiff.jacobian(u -> StabilizedNDEs.invariants(u, nothing, system, Val(2)), u)
    jac_analytic_2 = StabilizedNDEs.invariants_jacobian(u, nothing, system, Val(2))
    @test jac_analytic_2 ≈ jac_forward_diff_2 atol = 1e-6
end
