% Newton Method as a subroutine. Modified on Sep 24, 2004
% Use this progrqam in combination with an auxiliary program
% called TestNewt which is also on my website. In Test Newt
% supply all the parameter values. 
% Also, the functions 'fn' and 'fnprima' are on my website.


function Newt (ft,fnprima,Tol,itmax, A, B,p0)

format short e
% **************************INPUT DATA***********************************.

disp('   Input data    ');
disp('     A              B          TOL          itmax          po     ');
disp( [A,B,Tol, itmax,p0] );

%*******GRAPH OF THE FUNCTION FOR WHICH ROOTS ARE BEING FOUND*************

% Partition of the interval x used to obtain the graph of f(x).

step=(B-A)/100;
x1=A:step:B;
y1=feval(ft,x1);
plot(x1,y1)
xlabel('x'); ylabel('f(x)');
grid on

% *******************MATLAB APPROXIMATION OF ROOT************************

MATLABApprox=fzero(ft,p0);
format long e;
disp('MATLAB Approx.');
disp(MATLABApprox);

%*******CHECKING IF THERE IS A ROOT INSIDE THE INTERVAL SELECTED*********
format short e;
flag=0;  
if feval(ft,A)*feval(ft,B)>0
    disp('There is not root in the given interval');
   	disp(' Choose new values for the first two approximations');
    return
end
if feval(ft,A)==0 
    disp('The left end A is an exact root');
    fprintf('Exact root = %f\n',A); 
    flag=1;
elseif feval(ft,B)==0
    disp('The right end A is an exact root');
    fprintf('Exact root = %f\n',B);
    flag=1;
end

% *********************************MAIN  LOOP****************************
disp('         it        approx.      ErrBound       f(pn)     AbsErr  ');                 
it=0;
while it<itmax 
    it=it+1;
    pn = p0 - feval(ft,p0)/feval(fnprima,p0);
    ErrBound=abs(pn-p0);
    AbsErr=abs(MATLABApprox-pn);
    disp([it,pn,ErrBound,feval(ft,pn),AbsErr]);
    if feval(ft,pn)==0 
        disp('pn is an exact root');
        disp('Exact Root and Number of Iterations');
        format long e;
        disp([pn,it]);
        return
    end
% STOP CRITERIA (Comparing two consecutive approximations)
    if ErrBound < Tol
        disp('Approximated Root and Number of Iterations');
        format long e;
        disp([pn,it]);
        return
    else 
        p0=pn;
    end
end
disp('More iterations are needed');