subroutine try_newton_solve(n, x0, residual_fn, x_out, final_norm, success)
integer, intent(in) :: n
real(dp), intent(in) :: x0(n)
procedure(nonlinear_residual) :: residual_fn
real(dp), intent(out) :: x_out(n)
real(dp), intent(out) :: final_norm
logical, intent(out) :: success
integer :: iter, backtrack
real(dp) :: x(n), f(n), jac(n, n), dx(n), x_trial(n), f_trial(n), step_scale, fnorm, trial_norm
logical :: linear_ok, improved
x = x0
call residual_fn(x, f)
fnorm = residual_norm(f)
do iter = 1, nonlinear_max_iter
if (fnorm < nonlinear_tol) exit
call numerical_jacobian(n, x, f, residual_fn, jac)
call solve_small_linear_system(n, jac, -f, dx, linear_ok)
if (.not. linear_ok) exit
step_scale = 1.0d0
improved = .false.
do backtrack = 1, nonlinear_max_backtrack
x_trial = x + step_scale*dx
call residual_fn(x_trial, f_trial)
trial_norm = residual_norm(f_trial)
if (trial_norm < fnorm) then
x = x_trial
f = f_trial
fnorm = trial_norm
improved = .true.
exit
end if
step_scale = 0.5d0*step_scale
end do
if (.not. improved) exit
end do
x_out = x
final_norm = fnorm
success = fnorm < nonlinear_tol
end subroutine try_newton_solve