PRINT; -reset false; -user_print true
SOLUTION 1
USER_PRINT
# find x that gives f_x = 0 with Newton iteration.
# Here, f_x = x^0.5 - 2^0.5. Thus, x = 2.
# You can include the lines 9400-9880 in a PHREEQC input file and another USER_ block.
#   Define f_x in line 9410.
#   Set options in lines 9510-9550.

#   Give an initial value for x, and then, GOSUB 9500.

10 x = 0 # initial value for x.
20 GOSUB 9500
# 30 if abs(x - LB_x) <= precision then print 'WARNING: x is equal to lower bound, please check result...'
# 40 if abs(x - UB_x) <= precision then print 'WARNING: x is equal to upper bound, please check result...'
50 print 'Result = ', TRIM(STR$(x)), 'after', TRIM(STR$(it)), 'iterations.'
60 end

#  Define f_x =======================
9410 f_x = x^0.5 - 2^0.5

#  Set options ======================
9510 LB_x = 0             # estimated lower bound for x
9520 UB_x = 20            # estimated upper bound for x
9530 precision = 1e-10    # iterate until this error
9540 max_it = 25          # maximal iteration steps
9550 deriv_x = 1e-3       # dx = deriv_x * x, central differences over 2 dx.

# program lines =====================
9400 rem # f_x
# 9410    f_x = x^0.5 - 2^0.5
9420 return

9500 rem # Newton(f_x, x)
9600  if x <= LB_x then x = LB_x + 2 * deriv_x
9610  if x >= UB_x then x = UB_x - 2 * deriv_x
9620  GOSUB 9400
9630  it = 0
9640  while abs(f_x) > precision and it < max_it
9650    it = it + 1
9660    xn1 = x
9670    GOSUB 9400
9680    f_x1 = f_x
# 9690  print it, x, f_x
9700    GOSUB 9800
9710    a = xn1 - f_x1 / deriv_f_x
9720    # do some checks on a
9730    if a < LB_x then a = LB_x : if it = max_it then print eol$ + 'WARNING: x is at lower bound =', trim(str$(a)), '. Please decrease lower bound...' + eol$
9740    if a > UB_x then a = UB_x : if it = max_it then print eol$ + 'WARNING: x is at upper bound =', trim(str$(a)), '. Please increase upper bound...' + eol$
9750    x = a
9760    GOSUB 9400
9770  wend
9780  return

9800 rem # calculate the derivative of f_x
9810   if x = 0 then dx = 2 * deriv_x else dx = deriv_x * x
9820   x = x + dx
9830   GOSUB 9400
9840   f_xn1 = f_x
9850   x = x - 2 * dx
9860   GOSUB 9400
9870   deriv_f_x = (f_xn1 - f_x) / (2 * dx)
9880  return
-end
END