function [x,iter] = sqrtn(a,tol)
%SQRTN Square root of a scalar by Newton's method.
% X = SQRTN(A,TOL) computes the square root of the scalar
% A by Newton's method (also known as Heron's method).
% A is assumed to be >= 0.
% TOL is a convergence tolerance (default EPS).
% [X,ITER] = SQRTN(A,TOL) returns also the number of
% iterations ITER for convergence.
if nargin < 2, tol = eps; end
x = a;
iter = 0;
xdiff = inf;
fprintf(' k x_k rel. change\n')
while xdiff > tol
iter = iter + 1;
xold = x;
x = (x + a/x)/2;
xdiff = abs(x-xold)/abs(x);
fprintf('%2.0f: %20.16e %9.2e\n', iter, x, xdiff)
if iter > 50
error('Not converged after 50 iterations.')
end
end