Hi there,

I have a multivariate rootfinding problem, where one of the functions utilizes `QuantEcon.jl`

's `qnwnorm`

which generates the nodes for the Gauss-Hermite quadrature. However, using `NLsolve`

on one of the parameters that is an input to `qnwnorm`

throws an error, as `qnwnorm`

tries to transform the dual number to a `Float64`

, which fails.

A admittedly contrived MWE would be something like this:

```
using QuantEcon, NLsolve
function f!(res, x)
nodes, weights = qnwnorm(10,x[1],1)
res[1] = x[1] + x[2] - 3.5
res[2] = mean(nodes) - 1
return res
end
nlsolve(f!, [1., 1.], autodiff = :forward)
```

which throws

```
MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{typeof(f!),Float64},Float64,2})
Closest candidates are:
Float64(::Real, !Matched::RoundingMode) where T<:AbstractFloat at rounding.jl:194
Float64(::T<:Number) where T<:Number at boot.jl:718
Float64(!Matched::Int8) at float.jl:60
...
```

Is there a different solver/package in Julia that can handle these kinds of problems?

Thanks a lot!