using Plots
gr(size=(400, 300), fmt=:png);

using NLsolve

using ForwardDiff

using LinearAlgebra



a = Ref{Float64}(2.)

function f!(du, u)
    x, y = u
    du[1] = a[] / (1 + y^2) - x
    du[2] = a[] / (1 + x^2) - y
    nothing
end

a[] = 3.
r = nlsolve(f!, [1., 4.], autodiff=:forward)

r.zero



n = 10000

as = rand(n) * 6 # list of a
result = rand(n, 2) * 6 # use as initial u

@time for i in 1:n
    
    a[] = as[i]
    r = nlsolve(f!, result[i, :], autodiff=:forward)
    
    if r.f_converged
        result[i, :] = r.zero
    end
end

scatter(as, result[:, 2],
    markeralpha=0.2, markersize=2, markerstrokewidth=0,
    legend=:none, xlims=(0, 6), ylims=(0, 6))


scatter(as, result[:, 1], result[:, 2],
    xlims=(0, 6), ylims=(0, 6), zlims=(0, 6),
    markeralpha=0.2, markersize=2, markerstrokewidth=0,
    legend=:none)



negative(x) = x < 0
stable(J) = reduce(&, negative.(eigvals(J)))

a[] = 4.
@show r = nlsolve(f!, [3., 3.], autodiff=:forward)

@show J = ForwardDiff.jacobian(f!, similar(r.zero), r.zero)
@show eigvals(J)
@show stable(J)



stability = zeros(Int, n)
dummy = zeros(2)

@time for i in 1:n
    
    a[] = as[i]
    J = ForwardDiff.jacobian(f!, dummy, result[i, :])
    
    if stable(J)
        stability[i] = 1
    end
end

scatter(as, result[:, 2], color=stability,
    markeralpha=0.2, markersize=2, markerstrokewidth=0,
    legend=:none, xlims=(0, 6), ylims=(0, 6))