Use nonlinear solver in gamma_loop

[1-Bart-1]

According to Oriol (and me) using a nonlinear solver with jacobians should be more robust and need less iterations than the current method.

[1-Bart-1]

Using the unified interface for portability:
https://docs.sciml.ai/NonlinearSolve/stable/

[ufechner7]

Sound good in my ears!
According to experiance the Julia solvers outperform most other solutions. But we also need to look at the memory allocations.

[1-Bart-1]

Will use a pure julia solver with AutoFiniteDiff() first, as this is easiest to implement. The code is incompatible with ForwardDiff now, but finitediff should work well already.