%
% Call this function as shown below:
%
%    >> [x_bs i_bs x_nr i_nr]=rootfinding;
%    >> [x_bs' i_bs' x_nr' i_nr']
%    >> x_nr' / pi
%
function [x_bs i_bs x_nr i_nr] = rootfinding

for k=1:5 % k is the root (k-th) we are searching

    %
    % Binary search
    %
    % [We find the root approximately using BS with high tolorance]
    xmin = pi*(k-1) + pi*0.000001;
    xmax = pi*(k-1) - pi*0.000001 + pi/2;
    tol = 1e-2;

    fmax = func(xmax);
    fmin = func(xmin);

    i = 0;
    while (xmax-xmin)>tol
        i = i+1;
        xmid = (xmax+xmin)/2.0;
        fmid = func(xmid);

        if (fmax/fmid) > 0.0
            fmax = fmid;
            xmax = xmid;
        else
            fmin = fmid;
            xmin = xmid;
        end        
    end

    x_bs(k) = (xmax+xmin)/2.0;
    i_bs(k) = i; % The number of iteration of bin-search

    
    
    %
    % Newton-Raphson
    %
    % [Then we request a high precision solution starting with the
    %  solution found by the binary search]
    tol = 1e-6;
    x = x_bs(k);
    xpre = x + pi*0.001;

    f = func(x);
    fpre = func(xpre);

    i = 0;
    while (abs(x-xpre)>tol)
        i = i + 1;

        xnew = x  - f*(x-xpre)/(f-fpre);
        fnew = func(xnew);

        xpre = x;
        fpre = f;

        x = xnew;
        f = fnew;
    end

    x_nr(k) = x;
    i_nr(k) = i; % Number of iteration for N.R. search
end


%
% This part draw exp and tan seperately
% It takes care if divergence of tangent
% [Gets rid of +inf to -inf lines]
%
x = 0:pi/1000:pi*k;
xt = [x(diff(tan(x))>=0) x(diff(tan(x))<0)]';
tt = [tan(x(diff(tan(x))>=0)) inf(size(x(diff(tan(x))<0)))]';
o = [xt tt];
o = sortrows(o);
x = o(:,1);
t = o(:,2);
plot(x, exp(x/10), 'b', x, t,'g', x_nr, exp(x_nr/10),'ro')
axis([0 max(x) 0 1.1*max(tan(x_nr))])



%
% The function for which we are looking for the roots
%
function y = func(x)

y = exp(x/10) - tan(x);