Non-monotonic equations and `mcpsolver` from `NLsolve`

Dear all,

Working on some related topics, I’ve realized that the mcpsolver provided by the NLsolve package is struggling to find the solution(s) to non-monotonic functions.

As an example, say we are looking for the x’s which solve (x^3 + x^2 - 20x - 8)/6 - 3 = 0, subject to the constraint x >= 0. Graphically, this looks like

The three roots are complex-valued. However, one of them has a real part greater than 0, `x = 4.59056`, at which point the function evaluates at 0. Tackling the problem with the `mcpsolver`, we have
function f!(x,fvec) 
  fvec[1] = (x[1]^3 + x[1]^2 - 20x[1] - 8)/6 - 3
end

NLsolve.mcpsolve(f!, [0.], [Inf], [0.], reformulation = :smooth, autodiff = true) # returns -1.48382
NLsolve.mcpsolve(f!, [0.], [Inf], [0.], reformulation = :minmax, autodiff = true) # returns -1.0
NLsolve.mcpsolve(f!, [0.], [Inf], [0.], reformulation = :smooth, autodiff = true, method = :newton) # returns 6.25915e8 

First, the solver returns negative value, even though the constraint indicates that x needs to be positive. Second, it does not find the root at 4.6.

I was wondering whether this is because the root is actually complex-valued. And if not, whether there would be a way to find this root, without changing the initial condition.

Thanks a lot for the help.

Have you looked at Roots.jl?

I had not because I was first working on multidimensional problems. But for functions with a single variable, like this one, Roots works indeed nicely. Thanks!