Newton solver in Stan

The Jacobian computation question aside, I propose two sets of optional tuning parameters. Available signatures would then be:

algebra_solver_newton(system, y_guess, theta, x_r, x_i);

algebra_solver_newton(system, y_guess, theta, x_r, x_i, rel_tol, 
                      abs_tol, max_steps);

algebra_solver_newton(system, y_guess, theta, x_r, x_i, 
                      rel_tol, abs_tol, max_steps,
                      line_search, inexact_jacobian);

In my view, these options remain fairly manageable. The first PR will cover the first two options, then I can add the third signature in a second PR. Hopefully, we’ll get feedback from users on what they need.