%**************************************************************************
%*     Problem 4.2                                                        *
%*     filename: problem04_2.m                                               *
%*     Programed by Ryoichi Kawai for Computational Physics Course        *
%*     Revised on 02/08/2016.                                             *
%**************************************************************************

clear all;

z0=6;  % system parameter

f=@(z) z*cot(z) + sqrt(z0*z0-z*z);  % define function

K=ceil(z0/pi);  % number of roots
N=100; 

% Plot the function and visually find blackets.
x0=0;
x1=pi;
for k=1:K
    dx = (x1-x0)/N;
    
    for n=1:N-1
        x(n) = x0 + dx*n;
        fx(n) = f(x(n));
    end
    
    p=plot(x(1:N-1),fx(1:N-1));
    set(p,'Linewidth',2);
    hold on
    x0 = x1;
    x1 = min([z0,x0 + pi]);
    
end
q=plot([0,z0],[0,0]);
set(q,'Color','black');
axis([0 z0 -10 10]);
r=plot([pi,pi],[-10,10],'--',[2*pi,2*pi],[-10,10],'--');
set(r(1),'Color','black');
set(r(2),'Color','black');
hold off

% initial bracket
z1=0;
z2=pi;

for k=1:K   
    dz = (z2-z1)/N;  % small shit
    z = rootfinding(f,z1+dz,z2-dz,10,10^(-6));
    fprintf('root=%.8f\n',z);
    z1 = z2;
    z2 = min([z0,z1 + pi]);
end
