gistfile1.txt

using PseudoArcLengthContinuation, LinearAlgebra, Plots, PyPlot, Optim
const PALC = PseudoArcLengthContinuation
using ForwardDiff

function f(x, alpha)
    [alpha*x[1] + x[1]^3 ]
end
function Jf(x, alpha)
    Jf = zeros(1, 1)
    Jf[1] = alpha+3x[1]^2
    Jf
end

Fmit = f
D(f, x, p, dx) = ForwardDiff.derivative(t->f(x .+ t .* dx, p), 0.)
d1Fmit(x,p,dx1) = D((z, p0) -> Fmit(z, p0), x, p, dx1)
d2Fmit(x,p,dx1,dx2) = D((z, p0) -> d1Fmit(z, p0, dx1), x, p, dx2)
d3Fmit(x,p,dx1,dx2,dx3) = D((z, p0) -> d2Fmit(z, p0, dx1, dx2), x, p, dx3)
jet = (f, Jf, d2Fmit, d3Fmit)


opt_newton = PALC.NewtonPar(tol = 1e-8, verbose = true)

# options for continuation
alpha_min = -10.
alpha_max = 10.
alpha_start = 2.
opts_br = ContinuationPar(dsmin = 0.001, dsmax = 0.05, ds = -0.01, pMax = alpha_max, pMin = alpha_min,
	detectBifurcation = 2, nev = 30, plotEveryNsteps = 10,maxSteps = 100, precisionStability = 1e-6, nInversion = 4, dsminBisection = 1e-7, maxBisectionSteps = 25)

optnewton = NewtonPar(verbose = true)

sol_start, _, _ = newton( x -> f(x, alpha_start), [1.], optnewton)

br, _ = continuation(f, sol_start, alpha_start, opts_br, plot = true, printSolution = (x,p) -> x[1], plotSolution = (x, p;kwargs...) -> (plot!(x; ylabel="solution",label="",kwargs...)))


br2, _ = continuation(jet..., br, 1, opts_br, plot = true, printSolution = (x,p) -> x[1], plotSolution = (x, p;kwargs...) -> (plot!(x; ylabel="solution",label="",kwargs...)))

Plots.plot([br,br2],plotfold=false)