Improve MESA's linesearch

MESA uses a Newton-Raphson scheme with line search to solve for the zero point of a linearized system and find the solution to the stellar structure. This involves inverting a matrix (with the bcyclic solver with pivoting and equilibration) and iteratively updating the solution until errors fall below prescribed tolerances.

See Section 6.3 of https://arxiv.org/abs/1009.1622

It's a known issue that MESA's solver can struggle in certain parts of parameter space where the system is stiff.

@Debraheem created a(/another) reproducer in https://github.com/MESAHub/mesa/pull/753 where @Debraheem and @pmocz are investigating precisely why the solver fails.

So far, we've found that when things break down the Jacobian has high condition number (>10^15), but the bcyclic solver is still able to solve the system with errors of O(10^-11).

There are indications the problem is in the line search algorithm:

star/private/star_solver.f90: subroutine adjust_correction

The function being minimized is oscillatory and the line search is not able to take the right sized steps to find the minimum. At present, the line research uses some heuristics and hard-coded parameters.

So an open research question is to improve the robustness of the linesearch.