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

%===================================================
% minimum of effective potential via Newton-Raphson
%===================================================
clear all 
close all
clc
format long

% the potential
Phi        = @(r)(2./r.^2-5./r);
Phi_deriv  = @(r)(-4./r.^3 + 5./r.^2);
Phi_dderiv = @(r)(12./r.^4 - 10./r.^3);

% 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')

%---------------------------------------------------------
% find minimum of potential via Newton-Raphson iterations
%---------------------------------------------------------
% initial guess for the minimum
rmin = 0.60;

DPhi_rmin = -4./rmin.^3 + 5./rmin.^2
n = 0;    % we also count the number of iterations

while abs(Phi_deriv(rmin))>eps
   rmin       = rmin - Phi_deriv(rmin)./Phi_dderiv(rmin);
   n          = n + 1;   
end

% or use your own function
[rmin, n] = NewRaph(Phi_deriv, Phi_dderiv, 0.6);

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

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