% Script: RJB_Newton3.m
% implement Newton's method, with analysis of method in mind
% first give function and interval, look at plot
%      then use method 
% iterates and error estimates plotted after method used
%
% Input to prompts on screen:  
%         fun is the function to be zeroed, input as string 'fname'
%         fprime is derivative of function, input as string 'fname'
%         a and b are the left and right end points for plotting
%         epsilon is the tolerance in the interval size (nonzero)
%         test is 1 if zero is in interval
%         epsilon is tolerance for error between successive iterates
%         feps is tolerance for function value
%
  close all
  test = 0;
  fun = input('Input the function name in quotes (no .m): ');
  fprime = input('Input the derivative name in quotes (no .m): ');
  while test ~= 1 
    a = input('Left endpoint for plot: ' );
    b = input('Right endpoint for plot: ' );
    x = linspace(a,b); y = feval(fun,x); plot(x,y,'-',x,0*x,'--');
    xlabel('x'); ylabel('y');
    test = input('If zero is in interval, enter 1, otherwise 0: ');
end
  p = input('Input initial point for Newton method: ' );
  % froot = input('Input root for model problem error caclulation:');
  epsilon = input('input nonzero tolerance for root:  ');
  feps = input('input nonzero tolerance for function:  ');
  f = feval(fun,p);  fp = feval(fprime,p);
  nmax = 20; 
  xiter = zeros(1,nmax); fiter = zeros(1,nmax); 
  xiter(1) = p; fiter(1) = f;
  n = 1;  xdiff = 1; 
% model sequence for comparison
  max1 = 10; nvalues = [1:max1];
  model = 2.^-(2.^nvalues);
% start method
  while(n <= nmax & ( abs(xdiff) >= epsilon | abs(f) >= feps ) )
      n = n + 1;
      p = p - f/fp;
      f = feval(fun,p);
      fp = feval(fprime,p);
      xiter(n) = p; 
      fiter(n) = f;
      xdiff = p - xiter(n-1);
  end
% compute error estimate using computed root value
% as approximation to root (only n-1 error values)
  nerr = (1:n-1);
  err = abs(diff(xiter));
%
% start plots of bisection behavior
  figure(2);
  semilogy(nvalues,model,'--')
  hold on;
  semilogy(nerr,err(1:length(nerr)),'-*')
  legend('model seq','|p_n-p|',0)
  xlabel('n')
  ylabel('p_{n}')
  hold off;
  pause
%
% next plot
% this time, Interval_n+1 vs Interval_n
  figure(3);
  loglog(abs(err(1:n-2)),abs(err(2:n-1)),'-*',model(1:max1-1),model(2:max1),'--')
  xlabel('E_{n-1}')
  ylabel('E_{n}')
  legend('Error estimate','Model sequence',0)
%
% next plot, function value vs n
%  figure(4);
%  semilogy(fiter)
%  xlabel('n')
%  ylabel('f(p_{n})')
%  pause
%
% last plot, new function value vs old function value
%  figure(5);
%  loglog(fiter(1:n-1),fiter(2:n))
%  xlabel('f(p_n)')
%  ylabel('f(p_{n+1})')