% Unit 3 (MATLAB Course @ UAM)
% (c) Alexander Knebe

%========================================================
% root of effective gravitational potential by bisection
% (using an anonymous function)
%========================================================
clear all 
close all
clc
format long

% the potential
Phi = @(r)(2./r.^2-5./r);

% start with simply plotting the potential
r   = linspace(0.2,2,1000);

plot(r,Phi(r),r,0), xlabel('r'), ylabel('Phi'), title('effective gravitational potential')

%-----------------------------
% root finding via bi-section
%-----------------------------
% initial guess for the root
a = 0.2;
b = 2.0;
n = 0;    % we also count the number of iterations

r0 = (a+b)./2;
while abs(a-b)>1e-8
   n  = n + 1;
   if Phi(a)*Phi(r0) > 0
       a = r0;
   else
       b = r0;
   end   
   r0 = (a+b)/2;
end

% or use your own function
[r0, n] = bisec(Phi, 0.2, 2.0);

r0             % print the root
n              % print the number of iteration


% mark root
hold on
plot([r0 r0],[min(Phi(r)) 0],'r')
