% This code is an implementation of Newton's method

% Specify the function f(x)

% NOTE: If you don't like/understand the syntax of the next 6 lines,
% you can also remove them, and specify the function f(x) as an
% external function (as before, cf. Bisectional or Fixed Point methods. 

fprintf(1,'Input the function f(x) in terms of x\n');
s = input(' ','s');
f = inline(s,'x');
fprintf(1,'Input the derivative of function f(x) in terms of x\n');
s = input(' ','s');
fprime = inline(s,'x');

fprintf(1,'Enter tolerance)\n');
tol = input(' ');

fprintf(1,'Enter maximum number of iterations)\n');
nmax = input(' ');

fprintf(1,'Enter initial guess)\n');
p0 = input(' ');

pold = p0;
false = 1;
true = 0;
finish = false;
iteration = 0;
epsilon = 1.0e-16;

% Beginning of the main loop
%%%%%%%%%%%%%%%%%%%%%%%%%%%%

while (finish == false & iteration <= nmax );

   iteration = iteration + 1;

   if (abs(fprime(pold)) < epsilon) ;
     fprintf(1,'derivative of f too close to zero; program will stop!');
     finish = true;
   else
     
     p = pold - (f(pold)/fprime(pold));

     fprintf(1,'iteration number %3d ',iteration);
     fprintf(1,'p is %15.8e \n',p);

     if (abs(p-pold)/abs(p) < tol) | (f(p)==0.0) ;

        % finished !

        p_final = p;
        finish = true;

     else;

        % not converged !

        pold = p;

     end

   end

end

if (iteration < nmax);

  fprintf(1,'result after iteration number %3d ',iteration);
  fprintf(1,'is %15.8e \n',p_final);

else

  fprintf(1,'no convergence after iteration number %3d ',nmax);

end
