% Function Newton-Raphson
%
% x0 : the initial estimation of the root
% acc : digits of the required accuracy
% f  : the function, f(root) ~= 0 
% df :  the direvative of the function
% root : the root
% no_fun : number of function evaluations
% tot : total number of repetitions
%
function [root,no_fun,tot] = newton_raphson(x0, acc, maxiter, f, df)

fout = fopen('results_newton_raphson.dat', 'w+');
fprintf(fout, ' Iter \t || \t x(k) \t|| \t f(x(k)) \t || \t df(x(k)) \t || |x(k)-x(k-1)| \n');   

root = x0;
tot = 1;
x(1) = x0;
fx(1) = feval(f,x(1));
dfx(1) = feval(df, x(1));
no_fun = 1;
accur=0.5*10^(-acc);

k = 1;

fprintf(fout, ' %5g \t || %10.5f || %10.5f || %10.5f || --- \n',k,x(k),fx(k),dfx(k));

while ( k <= maxiter )

    k = k + 1;
    x(k) = x(k-1) - fx(k-1)/dfx(k-1);%epanaliptiko sxima N-R
    fx(k) = feval(f,x(k));
    dfx(k) = feval(df, x(k));
    no_fun = no_fun + 2;
    
    dif = abs(x(k)-x(k-1));
    
    fprintf(fout, ' %5g \t || %10.5f || %10.5f || %10.5f || %10.5f \n',k,x(k),fx(k),dfx(k),dif);
    
    if ( (dif <= accur) && (abs(fx(k)) <= accur) )%kritirio termatismou tou algorithmou
        
        root = x(k);
        tot = k;
        fprintf(fout, ' The root is calculated after %5g iterations and it is: %10.5f \n', tot,root);
        fprintf(fout, ' Number of function evaluations = %5g \n', no_fun);
        break;
        
    end
    
end

fclose('all')
% Plot the graphical description of Newton-Raphson
%graph_newt(x,fx,'f');

end