% Iteration Method - Steffensen's Method - Newton's Method

g = inline('sqrt(10/(x+4))');
f  = inline('x-sqrt(10/(x+4))');
fp = inline('1+sqrt(10/(x+4)^3)/2');
p0 = 1.5;

tol      = 10^(-5);
MAX_ITER = 25;
i = 0;

fprintf('\n');
fprintf('-----------------Iteration---------------------------------------------\n');
fprintf(' i            pn             g(p)\n');
fprintf('-----------------------------------------------------------------------\n');

dx=2*tol;
while( i < MAX_ITER && abs(dx)>=tol ) 
  fprintf('%2d  % 8.5e  % 8.5e \n', i, p0, g(p0));
  i = i+1;
  p = g(p0);
  dx = p - p0;
  if( abs(dx)<tol )
    fprintf('%2d   %8.5e  % 8.5e\n',i, p0, g(p0));
  else
    p0 = p;
  end
end
if( i>= MAX_ITER )
  fprintf('WARNING: Fixed Point iteration did not converge after %d steps.\n',i);
end

x0 = 1.5;
x  = x0;
i  = 0;
dx=2*tol;

fprintf('---------Newton------------------\n');
fprintf('  i              x           f(x)\n');
fprintf('---------------------------------\n');
while( i < MAX_ITER && abs(dx)>=tol )
  fprintf('%3d   %8.5e  % 8.5e\n',i, x, f(x));
  dx = f(x)/fp(x);
  i = i+1;
  x = x - dx;
  if( abs(dx)<tol )
    fprintf('%3d   %8.6f   % 12.6e\n',i, x, f(x));
  else
    p0 = p;
  end
end
if( i>= MAX_ITER )
 fprintf('WARNING: Newton iteration did not converge after %d steps.\n',i);
end

i = 0;
p0 = 1.5;

fprintf('\n');
fprintf('-----------------Steffenson--------------------------------------------\n');
fprintf(' i            pn            p1           p2');
fprintf('             p          dx \n');
fprintf('-----------------------------------------------------------------------\n');
while( i < MAX_ITER )
  p1 = g(p0);
  p2 = g(p1);
  p  = p0-(p1-p0)^2/(p2-2*p1+p0);
  dx = p - p2;
  fprintf('%2d  % 8.5e  % 8.5e  %8.5e  % 8.5e  % 8.5e\n', ...
	  i, p0, p1, p2, p, dx);
  if( abs(dx)<tol )
    fprintf('%2d  % 8.5e  % 8.5e  %8.5e  % 8.5e  % 8.5e\n', ...
	  i, p0, p1, p2, p, dx);
    return;
  end
  i = i+1;
  p0 = p;
end

%-----------------Iteration---------------------------------------------
% i            pn             g(p)
%-----------------------------------------------------------------------
% 0   1.50000e+000   1.34840e+000 
% 1   1.34840e+000   1.36738e+000 
% 2   1.36738e+000   1.36496e+000 
% 3   1.36496e+000   1.36526e+000 
% 4   1.36526e+000   1.36523e+000 
% 5   1.36523e+000   1.36523e+000 
% 6   1.36523e+000   1.36523e+000
%---------Newton------------------
%  i              x           f(x)
%---------------------------------
%  0   1.50000e+000   1.51600e-001
%  1   1.36495e+000  -3.11226e-004
%  2   1.36523e+000  -1.35587e-009
%  3   1.365230    0.000000e+000
%
%-----------------Steffenson--------------------------------------------
% i            pn            p1           p2             p          dx 
%-----------------------------------------------------------------------
% 0   1.50000e+000   1.34840e+000  1.36738e+000   1.36527e+000  -2.11115e-003
% 1   1.36527e+000   1.36523e+000  1.36523e+000   1.36523e+000  -5.69959e-007
% 1   1.36527e+000   1.36523e+000  1.36523e+000   1.36523e+000  -5.69959e-007