Hi everyone! I’m relatively new to Julia (1.7.2), coming from a Matlab and Python background. I really enjoy the language, but I also struggle sometimes in cases where in Matlab / numpy I would simply use vectorization and it would be relatively fast.

I’m playing around with a nonlinear equation, the dispersion relation for ocean waves, which I’m trying to solve for k:

\omega^2 = gk \cdot tanh(kd)

In the formulation of my problem, \omega and k are vectors of length N and d is a scalar. I am using NonlinearSolve and setup the problem such that constant parameters in the optimization are provided in the NamedTuple p:

```
using NonlinearSolve
function dispersion(k::AbstractVector{<:Real}, p::NamedTuple)
@. k * p.G * tanh(k * p.d) - p.ω^2
end
# Setup parameters
N = 1000
G, d, ω = 9.8, 30.0, 2π * collect(range(0.05, 1, N))
p = (G=G, d=d, ω=ω)
# Starting point: first order approximation
k0 = @. ω^2 / G / tanh(ω * sqrt(d / G))
prob = NonlinearProblem{false}(dispersion, k0, p);
solver = solve(prob, NewtonRaphson(), tol=1e-9)
```

Doing some benchmarks against my original python implementation using scipy.newton shows that the julia version above dramatically underperforms, especially as N increases:

I played around with the tolerances as well as the Optim.jl package (which gave much worse results with different solvers). From the scipy docs I know that newton defaults to the secant method. I have not found a method yet that is faster than the NonlinearSolve approach above, but I’m also quite new to Julia, so I don’t know the entire landscape of solver packages yet :).

Does anybody know what causes the poor performance, or what solver package I should consider instead? Any thoughts are much appreciated!