Newton solver to find the mode of the Gaussian (for INLA)

[theorashid]

Part of INLA roadmap https://github.com/pymc-devs/pymc-experimental/issues/340.

We want to do the Laplace approximation around $\mu$. We need to find it first.

Notes

• jax implementation
• GPy implementation (maybe somewhere in here)