function [RootEstimateRange,Y] = newton2(x1,f,df,imax,tol,x2,n) 
% This function estimates multiple roots for equation 'f' for values...
% between x1 and x2. N is the number of intervals between to determine...
% guess for roots. 
% 'x' using the Newton Raphson method.
% Input: x1 initial lower guess at the root 
% f function to be evaluated 
% df derivative of f with respect to x 
% imax iteration maximum 
% tol convergence criterion 
% x2 Upper limit of guess
% n is number of intervals for linspace(x1,x2,n)
% Output: estimation of the root position  
% Newton's algorithm 
% xnew = xold + dx 
% dx = -fold(xold)/dfold(xold) 
 
x = linspace(x1,x2,n);

% initialise parameters:
x_next = inf;
diff = norm(x-x_next);
dfcheck = feval(df,x);
 
% CHECK: verify that the derivative is not zero. 
count=0; 
while diff>tol && count<=imax
    if dfcheck == 0           
        error('Problem');     
    end 

    f_at_guess = feval(f,x);
    dfdy_atfirst = feval(df,x);
    x_next = x - f_at_guess./dfdy_atfirst;
    diff=norm(x-x_next);

    count=count+1;       
    % CHECK: Excessive iterations 
    x = x_next;
end 
RootEstimate = unique(round(x_next,7));
Y = feval(f,RootEstimate);


% Used to Display and output roots only within x1 and x2 limits. 
[r,c] = size(RootEstimate);
RootEstimateRange(r,c) = NaN; 
for i=1:1:c
    if RootEstimate(1,i)<=x2 && RootEstimate(1,i)>=x1
    fprintf('\nRoot Estimates (6 d.p.): %.6f\n',RootEstimate(1,i))
    RootEstimateRange(1,i) = RootEstimate(1,i);
    else
     RootEstimate(1,i) = NaN;
    RootEstimateRange(1,i) = RootEstimate(1,i);
    Y(1,i) = NaN;
    end
end

RootEstimateRange = RootEstimateRange(~isnan(RootEstimateRange)); % remove NaNs
Y = Y(~isnan(Y)); % remove NaNs

end