function [y,ai]=newt(aa,bb);
close all
%
% p = newt(x0) USES NEWTON'S METHOD TO APPROXIMATE
% THE ROOT OF THE FUNCTION f(x) NEAR x0.
% THE ANSWER IS ACCURATE TO WITHIN 1E-5.
%
[F_plot,dF_plot]=Fnewt([aa:(bb-aa)/100:bb]);

plot([aa:(bb-aa)/100:bb],F_plot);
hold on
plot([aa:(bb-aa)/100:bb],dF_plot);
legend('F','derivative')
drawnow

flag=input('Happy with interval?')
while flag==0
    aa=input('lower bound')
    bb=input('upper bound')
    close 
    [F_plot,dF_plot]=Fnewt([aa:(bb-aa)/100:bb]);
plot([aa:(bb-aa)/100:bb],F_plot);
hold on
drawnow
flag=input('Happy with interval?')

end
x0=input('Initial guess')

tol = 1e-10; N0 = 25;

y(1) = x0;

i = 2;
while i <= N0,

  [fm,dfm]=Fnewt(x0);
plot(x0,fm,'o')
drawnow
pause
  xm = x0 - fm/dfm; 
  y(i) = xm;
  if abs(xm-x0) < tol,
    p = xm; break;
  end
  i=i+1;
  x0 = xm;
end
y = y';
  p = xm;
  [fm,dfm]=Fnewt(x0);
  xm = x0 - fm/dfm; 

  
disp(' ')
disp(['The root is ',num2str(p,8),' with a possible error of ',num2str(abs(xm-x0))])
disp(' ')
if i> N0,
   disp(['Newton''s method did not meet the'])
   disp(['tolerance ',num2str(tol),' after ',num2str(N0),' iterations']),
   disp(' ')
end
% %-------- Aiken -----------
% for kk=1:length(y)-2
%    ai(kk) = y(kk) - (y(kk+1)-y(kk))^2/(y(kk+2)-2*y(kk+1)+y(kk));
% end
% ai(length(y)-1) = ai(length(y)-2);
% ai(length(y)) = ai(length(y)-2);
% ai = ai(:);
% 
% 
% err = [abs(  y - y(end)  )];  %err = [abs( ai-ai(end) )];
% x1 = err(1:end-1);  y1 = err(2:end);  
% plot(log(x1),log(y1)); hold on, plot(log(x1),log(y1),'r.'),hold off
% axis equal,grid
% 
% disp(' n         p_n             Aiken(p_n)         F(p_n)')       
% disp('----------------------------------------------------')
% for kk=1:length(y)
% disp(sprintf('%2.0f   %15.12f   %15.12f   %15.12f', kk, y(kk), ai(kk), Fnewt(y(kk))))
% end
% disp(' ')
% 
%--------------------------------
function [Fun,dFun] = Fnewt(x,h,ti,wi)

%Fun  = x.^6 - x - 1;  dFun = 6*x.^5 - 1;
Fun  = x.^3 + 4*x.^2 - 10; dFun = 3*x.^2 + 8*x;
%Fun  = sin(x);  dFun = cos(x);
%Fun = x - 0.2*sin(x) - 0.5;  dFun = 1 - 0.2*cos(x);
%Fun = cos(x)-x;   dFun = -sin(x)-1;
%Fun = exp(x)-2*x-1;   dFun = exp(x) - 2;
%Fun = exp(x) - x -1;  dFun = exp(x) -1;
%Fun = 16*x.^4-40*x.^3+5*x.^2+20*x+6;  dFun = 48*x.^3-120*x.^2+10*x+20;
%Fun = x - h*ti*exp(-x) - wi;  dFun = 1 + h*ti*exp(-x);
%Fun = x.^2+x+1;  dFun = 2*x + 1; %%%%%%%
%Fun = 230*x.^4 + 18*x.^3 + 9*x.^2 - 221*x -9; dFun= 920*x.^3+54*x.^2+18*x-221;