function x = newton_sqrt(a,atol)
% This function computes the square root x of the input a, assumed
% a non-negative number. It does so by applying Newton's method to the
% equation x^2 - a = 0, which is also known as the Babylonian method.
% The iteration is given by x_{n+1} = 0.5 * (x_n + a/x_n).
% It is known to converge globally from any positive initial guess,
% see https://math.stackexchange.com/questions/82682.
% Convergence is reached when abs(x^2 - a) <= atol.
% Assign some non-negative initial iterate and evaluate the initial residual.
x = a/2;
r = x^2 - a;
% Within the Newton loop, update the iterate and residual once per
% iteration until the residual tolerance is reached.
while (abs(r) > atol)
x = 0.5 * (x + a / x);
r = x^2 - a;
end