%**************************************************************************
%*     Problem 4.1                                                        *
%*     filename: problem04_1.m                                               *
%*                                                                        *
%*     Programed by Ryoichi Kawai for Computational Physics Course        *
%*     Revised on 02/08/2016.                                             *
%**************************************************************************
clear all;

% define function and its derivative
f=@(x) exp(-2*x)-2*exp(-x)+0.5;
df=@(x) 2*exp(-x)*(1-exp(-x));

% From the derivative, it is obvious that the potential is minimum
% at x=0.  Therfore, we plot the potential from xmin to xmax where
% xmax is a positive value and xmin negagive.

N=100;
xmax=2.0
xmin=-1.0;
dx=(xmax-xmin)/N;
x=linspace(xmin,xmax,N+1);
y=f(x);
plot(x,y,[xmin,xmax],[-0.5,-0.5])

% By direct inspection, the first turinig point is between x=-1 and x=0.
% The second turning point is between x=1 and x=2.


% initial bracket
x1=-1;
x2=0;

% tolerance
epsilon=1e-8;

for i=1:2
    
    % check the bracket
    f1=f(x1);
    f2=f(x2);
    if f1*f2 > 0
       error('Bracket is incorrect');
    end
 
    
   % bisection 10 iterations

   % reset iteration counter
   n=0;

   % Bisection method (10 iterations)
   xm = (x1+x2)/2;
   fm = f(xm);
   while n<10
      if f1*fm < 0  % root in the lower half
         x2=xm;
         f2=fm;
      else          % root in the upper half
         x1=xm;
         f1=fm;
      end
      xm = (x1+x2)/2;  % new mid point
      fm = f(xm);
      n=n+1;
   end

   fprintf('Bisection      = %.8f  (iteration= %i)\n',xm,n);

   % Newton-Raphson method
   x = xm;
   fx = fm;
   n = 0;
   while abs(fx)> epsilon
      dfx = df(x);
      x = x - fx/dfx;
      fx = f(x);
      n = n+1;
   end

   fprintf('Newton-Raphson = %.8f  (iteration= %i)\n',x,n);

   xc(i)=x;
   
   % set the bracket for the second turning point
   x1=1;
   x2=2;
end

fprintf('1st turning point = %9.6f\n',xc(1))
fprintf('2nd turning point = %9.6f\n',xc(2))
