"""

`freqamp` frequency perturbation amplitude
`n` number of perturbations 

Example
≡≡≡≡≡≡≡≡

julia> 
    J = 1       # only node 1 
    J = [1,2,3] # set of nodes 

"""
function stability_problem(p, J; freqamp = 1.0, t_final=120.0, tolerance = 0.1)

    # default ode problem
    u_sync = synchronous_state(p)
    ode_prob = ODEProblem(swing_dynamics!, u_sync, (0, t_final), p)

    # low discrepancy sequence for unit inverval
    N = nv(p.grid)
    m = length(J)
    P = primes(4*m^2+1)[1:2m]
    H = Halton.(P)

    # perturb intial condition 
    function _bstab_init(prob, i, repeat)
        u0 = u_sync[:]
        Z = getindex.(H, i) .* 2 .- 1
        u0[J]      .= Z[1:m] .* pi
        u0[J .+ N] .= Z[1+m:2m] .* freqamp
        prob.u0 .= u0
        prob
    end

    # check convergence to sync
    function _bstab_out(sol, i) 
        u = sol.u[end];
        Δu = maximum( abs.( u[1+N:end] .- u_sync[1+N:end] ) )
        S = (Δu < tolerance)
        return (S, false)
    end

    return prob = EnsembleProblem(ode_prob; 
        output_func = _bstab_out,
        prob_func = _bstab_init
    )    
end