% sqrt_fn.m 
%
% This functions computes the squares root of a number R numerically using
% Newton's method applied to 
%
% f(beta) = beta^2-R = 0
%
% or, analogously, applied to the problem
%
% argmin_beta {F(beta)=beta^3/3-R*beta}
%
% Note that F'(beta)=f(beta).
%
clear all, close all
R = input(' Enter a positive number R = ');
beta = R;
iterates = R;
iter_flag = 1;
while iter_flag
  % Newton iteration
  beta = 0.5*(beta + R/beta);
  % Story iteration history and check stopping criteria
  iterates = [iterates, beta];
  if beta^2-R < 1e-10
      iter_flag = 0;
  end
end
fprintf('The square root of R = %2.9f\n',beta);

% Plot iterates, function for intuition
beta_vec = linspace(0,R);
figure(1)
  plot(beta_vec,beta_vec.^2-R)
  hold on, plot(beta_vec,beta_vec.^3/3-R*beta_vec,'k')
  plot(iterates,zeros(size(iterates)),'r.')
  plot(beta_vec,zeros(size(beta_vec)),'g'), hold off
  legend('f(\beta)=\beta^2-R','F(\beta)=\beta^3/3-R*\beta','Newton iterates')