ENH: optimize._hessian: use _jacobian to compute accurate Hessian

[mdhaber]

Reference issue

Follow-up to gh-20630

What does this implement/fix?

Adds a private function to compute accurate Hessians using a finite difference approximation. Vectorized to compute Hessian at an arbitrary number of points in a single call.

Additional information

This is a pretty rough draft. We need to:

• Test (and fix, where applicable) attributes of return object (for this and _jacobian). It doesn't look like nfev is correct.
• Improve error estimate by considering error of _jacobian. (@steppi)
• Consider adding features to help user vectorize with respect to multiple parameters. (Probably a separate PR.)
• Consider using first Jacobian evaluation to estimate appropriate step size. This could reduce number of iterations required in subsequent evaluations. (Probably a separate PR.)

[mdhaber]

Not using callback right now.

Suggested change

	- ``-4`` : Iteration was terminated by `callback`.
	- ``1`` : The algorithm is proceeding normally (in `callback` only).