% compute square root of A using fixed point iteration
A  = 2;


% function defining the equation
f = @(x) x.^2 - A;
X = linspace(-2,2,200);
figure(1);
plot(X,f(X));
grid on;
xlabel('x'); ylabel('y');


% function defining the fixed point iteration
version = 0;
switch version
  case 0
    g  = @(x) A./x;
    dg = @(x) -A./x.^2;
  case 1
    % Richardson's iteration
    g  = @(x) x - (x.^2/A-1);
    dg = @(x) 1 - 2*x/A;
  case 2
    % Henron/babylonian/Newton's method
    g  = @(x) 0.5*(x+A./x);
    dg = @(x) 0.5*(1-A./x^2);
  case 3
    % Chebyshev's method
    g  = @(x) 0.375*x + 0.75*A./x    - 0.125*A^2./x.^3;
    dg = @(x) 0.375   - 0.75*A./x.^2 + 0.375*A^2./x.^4;
    d2g= @(x)           1.5 *A./x.^3 - 1.5  *A^2./x.^5;
end

figure(1);
X = linspace(1,2,200);
plot(X,X,X,g(X));
grid on;
xlabel('x'); ylabel('y');
legend('x','g(x)');

N    = 10;
tol  = 1.0E-08;
iter = 0;
x0   = 3;
r    = sqrt(A);
e0   = abs(x0-r);

while iter < N,
    iter = iter + 1;
    x1 = g(x0);
    e1 = abs(x1-r);
    
    fprintf('Iter: %3d, x: %15.10e, e1/e0: %e\n', iter, x1, e1/e0);

    if abs(x1-x0)/max(1.0,abs(x0)) < tol
        % success, converged!
        fprintf('\nConvergence!\n');        
        fprintf('Iter: %3d, x: %e, g'': %e\n', iter, x1, dg(x1));        
        break;
    end
    
    x0 = x1;
    e0 = e1;
end