% newtonraphson
function xr = newtonraphson(fun,dfun,x0)
% Demonstration of Newton-Raphson Method
%
% The input function should be of the form f(x)
% dfun is the derivative of the function f'(x)
%      
% Then call newtonraphson to find its roots
% xr = newtonraphson(@fun,@dfun,x0)
%
% For Example 4.8 use:
% newtonraphson(@(x) sin(5.*x)+cos(2.*x), @(x) 5.*cos(5.*x)-2.*sin(2.*x),0.5) 
%


close all
clc

pmin = -1; % x limits of the plot. Change for different functions if needed.
pmax = 1;

% First plot f(x)to see the behaviour
figure
hold on
fplot(@(x) fun(x),[pmin,pmax],'LineWidth',2)
fplot(0,[pmin,pmax],'-k')
ylim([-1.5,2]) % y limits of the plot. Change for different functions if needed.
xlabel('x')
ylabel('F(x)')
title('Newton-Raphson Demo')
grid on

% Initialize variables
es = 0.0005;
maxit = 100;

% prepare the output header
fprintf('\n  #         xi          er\n');

count = 0;  % iteration counter
er = 1;     % initial error value to enter loop
xi = x0;     % initialize xi and xold

fprintf('%3g %10.6g %10.4g\n',count,xi,er);

erplot = zeros(1,5);

% main loop
while (abs(er) > es && count<maxit)
    count = count +1;
    xold = xi;  % keep the previous estimate
    fxold = fun(xold); % will be used for plotting
    dfold = dfun(xold); % derivate at xold
    
    xi = xold - fxold/dfold;
  
    plot([xold,xold,xi],[0,fxold,0],'--','Color',[0.8500, 0.3250, 0.0980]);  % allows us to visualize the approximations
       
    er = (xi-xold)/xi;  % relative error
    erplot(count) = er;
   
    fprintf('%3g %10.6g %10.4g\n',count,xi,er);
%    pause  % Wait for key press
end

if count >= maxit
    disp(' ')
    disp(['Failed to converge in ',num2str(maxit),' iterations'])
    disp(' ')
    disp('--------------------------------------------')
    disp(' ')
    xr = NaN;
else
    disp(' ')
    disp(['Root has been found at x = ',num2str(xi)])
    disp(' ')
    disp('--------------------------------------------')
    disp(' ')
    xr = xi;
end

figure
hold on
fplot(0,[0,count],'-k')
plot(erplot(1:count),'LineWidth',2)
title('Relative Error')
xlabel('Iteration Number')
ylabel('Relative Error (er)')
hold off

