Amazing how easy this was!
Apologies for the rusty documentation (in progress).
"""
input: x: x[1] and x[2] amplitude of cos and sin mode of single harmonic anzats
input: p: angular frequency of the driving excitation
ouput: mismatchvec: residual vector of equations determining x[1] and x[2]
"""
function duffing_hb_singlemode(x,p)
A, B = x
omd = p
m=1; freq = .5; om0=2*π*freq; ga=0.1; stiffnlin=10.; F0 = 1.;
mismathvec = zeros(eltype(x),2)
mismathvec[1] = m*(om0^2-omd^2)*A + ga*omd*B + stiffnlin*.75*(A^3 + A*B^2)
mismathvec[2] = -ga*omd*A + m*(om0^2-omd^2)*B + stiffnlin*.75*(B^3 + A^2*B) - F0
return mismathvec
end
prob = BifurcationProblem(duffing_hb_singlemode, [0.,0.], 0., 1; record_from_solution = (x,p; k...) -> norm(x))
br = continuation(prob, PALC(), ContinuationPar(p_min = 0., p_max = 10.))
plot(br, xlabel="driving frequency", ylabel="single mode amplidude")
