noise with warmup

[oameye]

@variables t xa(t) xs(t);
@variables ω ωa ωs γ δ J λ α;

equations = [
    d(d(xs, t), t) +
    (ωs^2 - λ * cos(2 * ω * t)) * xs +
    γ * d(xs, t) +
    (α / 2) * (
        (2 * δ^2 + J^2) * xs^3 + 3 * J^2 * xa^2 * xs + J * δ * xa^3 - 3 * J * δ * xa * xs^2
    ),
    d(d(xa, t), t) +
    (ωa^2 - λ * cos(2 * ω * t)) * xa +
    γ * d(xa, t) +
    (α / 2) * (
        (2 * δ^2 + J^2) * xa^3 + 3 * J^2 * xs^2 * xa - J * δ * xs^3 + 3 * J * δ * xs * xa^2
    ),
]

system = DifferentialEquation(equations, [xs, xa])
add_harmonic!(system, xs, ω) # xs will rotate at ω
add_harmonic!(system, xa, ω) # xa will rotate at ω
harmonic_normal = get_harmonic_equations(system);

res = 200
δn(δ, J) = δ / √(δ^2 + J^2);
Jn(δ, J) = J / √(δ^2 + J^2);
δb = 0.003;
Jb = 0.005;
fixed_ωλ =
    Dict(
        ωa => √(1.0^2 + √(δb^2 + Jb^2)),
        ωs => √(1.0^2 - √(δb^2 + Jb^2)),
        γ => 0.007,
        α => 1.0,
        δ => δn(δb, Jb),
        J => Jn(δb, Jb),
    )
varied_ωλ = (ω => range(0.987, 1.013, res), λ => range(1e-6, 0.04, res))

result_ωλ_normal = get_steady_states(
    harmonic_normal,
    varied_ωλ,
    fixed_ωλ;
    show_progress = true,
);

res = 200
δn(δ, J) = δ / √(δ^2 + J^2);
Jn(δ, J) = J / √(δ^2 + J^2);
δb = 0.003;
Jb = 0.005;
fixed_ωλ =
    Dict(
        ωa => √(1.0^2 + √(δb^2 + Jb^2)),
        ωs => √(1.0^2 - √(δb^2 + Jb^2)),
        γ => 0.007,
        α => 1.0,
        δ => δn(δb, Jb),
        J => Jn(δb, Jb),
    )
varied_ωλ = (ω => range(0.987, 1.013, res), λ => range(1e-6, 0.04, res))

method = WarmUp(;start_parameters = fill(0.001+0.01*im, 8))
result_ωλ_normal = get_steady_states(
    harmonic_normal,
    method,
    varied_ωλ,
    fixed_ωλ;
    show_progress = true,
);