function area=numnewtons(option)
%NUMNEWTON  An animation showing the Newton's method algorithm.
%USE:   Type  numnewton  at the MATLAB prompt

%Peter Dunn
%21 July 2000

if nargin==0,

h = figure('Name','Numerical Newton''s Method', ...
			  'tag','tagNewton','NumberTitle','off');
   ConvData = { '0.0001', '0.00001', '10' } ;
   
   set( h, 'UserData', ConvData );
   

   enterframe = uicontrol('Style','Frame','Units','normalized',...
            'Position',[0.02 0.02 0.79, 0.30]);

   funtext = uicontrol('Style','text','Units','normalized',...
             'Position',[0.03 0.23 0.15,0.07],'String','Function:',...
             'FontSize',14,...
             'HorizontalAlignment','center');
   fun = uicontrol('Style','edit','Units','normalized',...
             'Position',[0.18 0.23 0.52 0.07],...
             'FontSize',14,...
             'ToolTip','The function whose zeros we find',...
             'String','x.^2','tag','tagFun' );

   drvtext = uicontrol('Style','text','Units','normalized',...
             'Position',[0.03 0.13 0.15 0.07],'String','Deriv:',...
             'FontSize',14,...
             'HorizontalAlignment','center');
   drv = uicontrol('Style','edit','Units','normalized',...
         'Position',[0.18 0.13 0.52 0.07],...
         'FontSize',14,...
         'ToolTip','The derivative of the function',...
         'String','2*x','tag','tagDrv' );

   lowertext = uicontrol('Style','text','Units','normalized',...
         'FontSize',14,...
         'Position',[0.03 0.03 0.15 0.07],'String','Min x:' );
   lower = uicontrol('Style','edit','Units','normalized',...
         'FontSize',14,...
         'ToolTip','Minimum x-value to display',...
         'Position',[0.18 0.03 0.08 0.07],'String','-0.4','tag','tagLower');

   uppertext = uicontrol('Style','text','Units','normalized',...
         'FontSize',14,...
         'Position',[0.26 0.03 0.14 0.07],'String','Max x:');
   upper = uicontrol('Style','edit','Units','normalized',...
         'FontSize',14,...
         'ToolTip','Maximum x-value to display',...
         'Position',[0.39 0.03 0.08 0.07],'String','0.6','tag','tagUpper');

   evalattext = uicontrol('Style','text','Units','normalized',...
         'FontSize',14,...
         'Position',[0.48 0.03 0.13 0.07],'String','Guess:');
   evalat = uicontrol('Style','edit','Units','normalized',...
         'HorizontalAlignment','left',...
         'FontSize',14,...
         'ToolTip','The starting value',...
         'Position',[0.61 0.03 0.09 0.07],'String','0.5','tag','tagEvalat');

   speedtextup = uicontrol('Style','text','Units','normalized',...
         'FontSize',14,...
         'Position',[0.7 0.27 0.1 0.04],'String',{'Slow'});
   speedtextlo = uicontrol('Style','text','Units','normalized',...
         'FontSize',14,...
         'Position',[0.7 0.03 0.1 0.04],'String',{'Fast'});
   speed = uicontrol('Style','slider','Units','normalized',...
         'FontSize',14,...
         'ToolTip','Changes animation speed',...
         'Position',[0.73 0.08 0.05 0.18],'Min',0,'Max',2,...
         'SliderStep',[0.1 0.25],  'Value',1,'tag','tagSpeed');

   optbtn = uicontrol('Style','pushbutton','Units','normalized',...
             'Position',[0.81 0.24 0.16 0.0733],...
             'FontSize',14,...
             'ToolTip','Set convergence options',...
             'String','OPTIONS','Callback','numnewton(50); ');
   gobtn = uicontrol('Style','pushbutton','Units','normalized',...
             'Position',[0.81 0.1666 0.16 0.0733],...
             'FontSize',14,...
             'ToolTip','Start Animation',...
             'String','GO','Callback','numnewton(4); ');
   pausebtn = uicontrol('Style','pushbutton','Units','normalized',...
             'Position',[0.81 0.0933 0.16 0.0733],...
             'FontSize',14,...
             'ToolTip','Pauses Animation',...
             'String','PAUSE','Callback','numnewton(5);','tag','PauseButton');
   quitbtn = uicontrol('Style','pushbutton','Units','normalized',...
             'Position',[0.81 0.02 0.16 0.0733],...
             'FontSize',14,...
             'ToolTip','Closes the program',...
             'String','QUIT','Callback','delete(gcf);');

   uimenu('Label','DEMONSTRATIONS','tag','tagDemoMenu');
   uimenu(findobj(gcf,'tag','tagDemoMenu'), ...
          'Label','Cycling', ...
          'Callback','numnewton(10);');
   uimenu(findobj(gcf,'tag','tagDemoMenu'), ...
          'Label','Cycling (2)', ...
          'Callback','numnewton(11);');
   uimenu(findobj(gcf,'tag','tagDemoMenu'), ...
          'Label','No root exists', ...
          'Callback','numnewton(12);');
   uimenu(findobj(gcf,'tag','tagDemoMenu'), ...
          'Label','No root exists (2)', ...
          'Callback','numnewton(13);');
   uimenu(findobj(gcf,'tag','tagDemoMenu'), ...
          'Label','Divergence', ...
          'Callback','numnewton(14);');
   uimenu(findobj(gcf,'tag','tagDemoMenu'), ...
          'Label','Divergence (2)', ...
          'Callback','numnewton(15);');
   uimenu(findobj(gcf,'tag','tagDemoMenu'), ...
          'Label','Bad starting value', ...
          'Callback','numnewton(16);');
   uimenu(findobj(gcf,'tag','tagDemoMenu'), ...
          'Label','Zero derivative', ...
          'Callback','numnewton(17);');
   uimenu(findobj(gcf,'tag','tagDemoMenu'), ...
          'Label','Zero derivative (2)', ...
          'Callback','numnewton(18);');
   uimenu(findobj(gcf,'tag','tagDemoMenu'), ...
          'Label','Convergence to wrong root', ...
          'Callback','numnewton(19);');
   uimenu(findobj(gcf,'tag','tagDemoMenu'), ...
          'Label','Convergence to wrong root (2)', ...
          'Callback','numnewton(20);');
   uimenu(findobj(gcf,'tag','tagDemoMenu'), ...
          'Label','Infinite derivative', ...
          'Callback','numnewton(21);');

   uimenu('Label','OPTIONS','tag','tagOptionsMenu');
   uimenu(findobj(gcf,'tag','tagOptionsMenu'),...
          'Label','Set convergence options',...
          'Callback','numnewton(50);');
   uimenu(findobj(gcf,'tag','tagOptionsMenu'),...
          'Label','Reset default convergence options',...
          'Callback','numnewton(51);');

	axes('Position',[0.10 0.38 0.85, 0.52]);
	axis( [ -0.4 0.6 0 1 ] )
   
   set( gca,'FontSize',16);
   
   zoom;

   %numnewton(4);

elseif nargin==1,

   if option == 3, %NUMBER REGIONS

      xhi = str2num(get( findobj(gcf,'tag','tagUpper'),'String'));
      xlo = str2num(get( findobj(gcf,'tag','tagLower'),'String'));
      xval = str2num(get( findobj(gcf,'tag','tagEvalat'),'String'));

      if (xhi-xlo) <= 0 ,
         watchoff;
         herr=errordlg('Upper limit must be greater than lower limit.',...
                       'Error!');
         set(herr,'WindowStyle','modal');
         return;
      end;

      if ~ ( ( xval>=xlo ) & ( xval<=xhi ) ),
         watchoff;
         herr=errordlg([ 'Must be evaluating derivative between ',...
                         'upper and lower x limits.'], 'Error!');
         set(herr,'WindowStyle','modal');
         return;
      end;

   elseif option == 4, %GO pressed

      watchon;

      delete(gca);

      numnewton(3);

      dx = 1.0e-12;
      f = get( findobj(gcf,'tag','tagFun'),'String');
      d = get( findobj(gcf,'tag','tagDrv'),'String');

      %Check if a derivative is given:
      noDrv=0;
      if ( isempty(deblank(d)) ),
         noDrv=1;
      end;
      
      evalat = str2num(get( findobj(gcf,'tag','tagEvalat'), 'String'));

      haxes = axes('Position',[0.10 0.38 0.85, 0.52]);

      xu = str2num(get( findobj(gcf,'tag','tagUpper'),'String'));
      xl = str2num(get( findobj(gcf,'tag','tagLower'),'String'));

      x = linspace(xl,xu,1000);

      eval(['ffn=',f,';'],'ffn=NaN;');

      if any(isnan(ffn)),
         watchoff;
         herr=errordlg(...
              ['Couldn''t evaluate the function!  Check the function ',...
              'and the evaluation limits.'],...
              'Error evaluating function');
         set(herr,'WindowStyle','modal');
         return;
      end;

%     Plot function
      hold on;
      h = plot( x, ffn, 'Linewidth',3 );

      ax = axis;
      ax(3) = min([0, ax(3)]);
      axis(ax);

      %Plot y=0
      plot([xl, xu],[0 0],'k-');

%     Evaluate function at x
      x = evalat; 
      eval( [ 'fx = ',f,';'], 'fx=NaN;');

      %plot(x, fx, 's','MarkerFaceColor','r', ...
      %                'MarkerEdgeColor','r',...
      %                'MarkerSize',10);
      plot(x, 0, 's','MarkerFaceColor','g', ...
                      'MarkerEdgeColor','k',...
                      'MarkerSize',10);

%     Display
      disp(' ');
      if ( noDrv==1 ),
      disp('                                     Approximate');
      end;
      disp('        X-Value        Function       Derivative');
      disp('        ~~~~~~~        ~~~~~~~~       ~~~~~~~~~~');

      fprintf('%15.5f %15.7f',x,fx);

      abserr = 1;
      fx = 1;
      errflag = 0;

      %Get convergence details
      ConvData = get( findobj('tag','tagNewton'), 'UserData');
      deltax = str2num( ConvData{1} );
      deltaf = str2num( ConvData{2} );
      maxnumits = str2num( ConvData{3} );

      if ( noDrv==1 ), 
         x = evalat - dx/2;
      end;
      numits = 0;
      while ( ( ( abs(fx)>deltaf ) & ( abs(abserr)>deltax ) ) & ...
           ( numits<maxnumits ) & ( errflag==0) )

         %Change colours of old lines
         if exist('plotnow','var'),
            set( plotnow, 'MarkerFaceColor',[0.4 0.4 0.4],...
                          'MarkerEdgeColor',[0.4 0.4 0.4],...
                          'MarkerSize',8);
            set( plotline, 'Color',[0.4 0.4 0.4]);
            set( plotdash, 'Color',[0.4 0.4 0.4]);
         end;
         
         % Indicate the guess
         eval( [ 'fx = ',f,';'], 'fx=NaN;');
         plotdash = plot( [x, x],[0,fx],'r:', 'LineWidth',2);

         %Determine tangent
x1=x;
         if ( noDrv==1), 
            %Numerically approximate the slope, m
            x = x+dx; 
x2=x;
            eval( [ 'fdx = ',f,';'], 'fdx=NaN; errflag=1;');
            m = ( fdx - fx ) / dx;

         else
            %Use the given derivative to get the slope, m

            ftmp = f;
            f = d;
            eval( [ 'm = ',f,';'], 'd=Nan; errflag=1;');
            d = f;
            f = ftmp;

         end;
         fprintf('%15.5f \n',m);

         %Check for evaluation problems
         if ( ( abs(m) < 1.0e-06 ) | ...
              ( isinf( abs(1/m) ) ) ), %Problems! Zero derivative!
            errflag = 1;
            %m = 0;
         end;
         if isinf(m),   %Infinite derivative
            errflag = 1;
         end;

         if ( errflag == 0 ) , 
         
            c = fx - m*x;
            nextzero = -c/m;
%[m x1 x2 fdx fx fdx-fx dx, nextzero]
            storex = x;
            x = nextzero; 
            eval( [ 'nextf = ',f,';'], 'nextf=NaN');
            x = storex;

            % Show the Tangent
            xvec = linspace(x, nextzero,2);
            tangentline = m*xvec+c;

            for i=2:length(xvec)

               plotline = plot([xvec(i-1),xvec(i)],...
                 [tangentline(i-1),tangentline(i)],'r-', ...
                 'LineWidth',2);
                 drawnow;

            end;

            %Show the next guess
            plotnow = plot(nextzero,0,'go','MarkerFaceColor','r',...
                           'MarkerEdgeColor','r',...
                           'MarkerSize',10,...
                           'tag','NextZero');

            lena = length(num2str(nextzero));
            lenb = length(num2str(nextf));
            lenm = length(num2str(m));
        
            fprintf('%15.5f %15.7f',nextzero,nextf);


            %Update for next iteration
            abserr = x - nextzero;
            x = nextzero;
            numits = numits+1;
            set(gca,'FontSize',16);

            %Check if speed has changed
            pausetime = get( findobj(gcf,'tag','tagSpeed'),'value');
            pause(pausetime);

         else

            watchoff;
            herr = errordlg( { 'Error in evaluation!  ',...
                               ' ',...
                               'Possible causes of the problem:',...
                               '* infinite values; ', ...
                               '* undefined function values; ',...
                               '* undefined derivative values; ',...
                               '* a zero derivative;',...
                               '* incorrect MATLAB function.',...
                               ' ',...
                               'Please check the inputs.'  },...     
                             'Error!' );             
            set( herr,'WindowStyle','modal');

         end;
 
         conv = 0;
         if ( abs(fx) <= deltaf ),
            disp(' ');
            disp('Convergence detected:  Function value close to zero.');
            conv = 1;
         end;
         if ( abs(abserr) <= deltax ),
            disp(' ');
            disp('Convergence detected:  x value not changing.');
            conv = 1;
         end;

      end;

      fprintf('\n');

      if (conv==0) & ( numits == maxnumits ),
         watchoff;

         errtext = ['No convergence after ',num2str(maxnumits),...
                    ' iterations!  Try a better starting point.'];
         disp(' ');
         disp(errtext)

         herr=errordlg(errtext, 'No Convergence!');
         set(herr,'WindowStyle','modal');

         return;
      end;
 
      hold on;

      %Font Size
      set(gca,'FontSize',16);

      drawnow;
         
      watchoff;
      hold off;

   elseif option==5, %Pause pressed

      %Test to see if watch is set to on
      pointertype = get( gcf, 'Pointer');
      if strcmp(pointertype,'watch') ,
         set( findobj(gcf,'tag','PauseButton'),'String','RESUME');
         set( findobj(gcf,'tag','PauseButton'),'Callback','numnewton(6);');
         set( findobj(gcf,'tag','PauseButton'),'ToolTip','Resumes animation');
         watchoff;
         waitfor( findobj(gcf,'tag','PauseButton'),'String','PAUSE');

         watchon;
      end;

   elseif option==6, %Resume pressed

      set( findobj(gcf,'tag','PauseButton'),'String','PAUSE');
      set( findobj(gcf,'tag','PauseButton'),'Callback','numnewton(5);');
      set( findobj(gcf,'tag','PauseButton'),'ToolTip','Pauses animation');

   elseif option==10, %Cycling 

      numnewton(51);

      set( findobj(gcf,'tag','tagUpper'), 'String','2');
      set( findobj(gcf,'tag','tagLower'), 'String','-3.5');
      set( findobj(gcf,'tag','tagFun'), 'String','x.^3-x-3');
      set( findobj(gcf,'tag','tagDrv'), 'String','3*x.^2-1');
      set( findobj(gcf,'tag','tagEvalat'), 'String','0');

      numnewton(4);

   elseif option==11, %Cycling (2)

      numnewton(51);

      set( findobj(gcf,'tag','tagUpper'), 'String','2');
      set( findobj(gcf,'tag','tagLower'), 'String','-2');
      set( findobj(gcf,'tag','tagFun'), 'String','sqrt(  abs(x) )');
      set( findobj(gcf,'tag','tagDrv'), ...
              'String','sign(x) ./ ( 2*sqrt(abs(x)) )');
      set( findobj(gcf,'tag','tagEvalat'), 'String','1');

      numnewton(4);

   elseif option==12, %No Root Exists

      numnewton(51);

      set( findobj(gcf,'tag','tagUpper'), 'String','6');
      set( findobj(gcf,'tag','tagLower'), 'String','0');
      set( findobj(gcf,'tag','tagFun'), 'String','x.^2 -4*x + 4.5');
      set( findobj(gcf,'tag','tagDrv'), 'String','2*x - 4');
      set( findobj(gcf,'tag','tagEvalat'), 'String','0.5');

      numnewton(4);

   elseif option==13, %No Root Exists (2)

      numnewton(51);

      set( findobj(gcf,'tag','tagUpper'), 'String','14');
      set( findobj(gcf,'tag','tagLower'), 'String','2');
      set( findobj(gcf,'tag','tagFun'), 'String', 'exp( -x )');
      set( findobj(gcf,'tag','tagDrv'), 'String','-exp( -x )');
      set( findobj(gcf,'tag','tagEvalat'), 'String','4');

      numnewton(4);

   elseif option==14,  %Diverging

      numnewton(51);

      set( findobj(gcf,'tag','tagUpper'), 'String','2');
      set( findobj(gcf,'tag','tagLower'), 'String','-4');
      set( findobj(gcf,'tag','tagFun'), 'String','atan(x)');
      set( findobj(gcf,'tag','tagDrv'), 'String','1 ./ (1+x.^2)');
      set( findobj(gcf,'tag','tagEvalat'), 'String','1.4');

      numnewton(4);

   elseif option==15,  %Diverging

      numnewton(51);

      set( findobj(gcf,'tag','tagUpper'), 'String','2');
      set( findobj(gcf,'tag','tagLower'), 'String','-3');
      set( findobj(gcf,'tag','tagFun'), 'String','erf(x)');
      set( findobj(gcf,'tag','tagDrv'), 'String','');
      set( findobj(gcf,'tag','tagEvalat'), 'String','1');

      numnewton(4);

   elseif option==16, %Bad starting value, all over the place

      numnewton(51);

      set( findobj(gcf,'tag','tagUpper'), 'String','2');
      set( findobj(gcf,'tag','tagLower'), 'String','-3');
      set( findobj(gcf,'tag','tagFun'), 'String','x.^3-x-1');
      set( findobj(gcf,'tag','tagDrv'), 'String','3*x.^2-1');
      set( findobj(gcf,'tag','tagEvalat'), 'String','-0.5');

      numnewton(4);

   elseif option==17, %Zero derivative

      numnewton(51);

      set( findobj(gcf,'tag','tagUpper'), 'String','5');
      set( findobj(gcf,'tag','tagLower'), 'String','-4');
      set( findobj(gcf,'tag','tagFun'), 'String','sin(x)');
      set( findobj(gcf,'tag','tagDrv'), 'String','cos(x)');
      set( findobj(gcf,'tag','tagEvalat'),...
            'String','4.36918237257240');

      numnewton(4);

   elseif option==18, %Zero derivative (2)

      numnewton(51);

      set( findobj(gcf,'tag','tagUpper'), 'String','2');
      set( findobj(gcf,'tag','tagLower'), 'String','-0.5');
      set( findobj(gcf,'tag','tagFun'), 'String', ...
           'x.^3-3*x.^2+3*x');
      set( findobj(gcf,'tag','tagDrv'),...
           'String','3*x.^2-6*x+3');
      set( findobj(gcf,'tag','tagEvalat'), 'String','1.7937');

      numnewton(4);

   elseif option==19, %Convergence to `wrong' root

      numnewton(51);

      set( findobj(gcf,'tag','tagUpper'), 'String','pi/2');
      set( findobj(gcf,'tag','tagLower'), 'String','-5*pi/2');
      set( findobj(gcf,'tag','tagFun'), 'String','sin(x)');
      set( findobj(gcf,'tag','tagDrv'), 'String','cos(x)');
      set( findobj(gcf,'tag','tagEvalat'), 'String','1.25');

      numnewton(4);

   elseif option==20, %Convergence to `wrong' root (2)

      numnewton(51);

      set( findobj(gcf,'tag','tagUpper'), 'String','1.5');
      set( findobj(gcf,'tag','tagLower'), 'String','-1');
      set( findobj(gcf,'tag','tagFun'), 'String','humps(x)');
      set( findobj(gcf,'tag','tagDrv'), 'String','');
      set( findobj(gcf,'tag','tagEvalat'), 'String',-0.55');

      numnewton(4);

   elseif option==21, %Infinite derivative

      numnewton(51);
      s = warning;
      warning off;

      set( findobj(gcf,'tag','tagUpper'), 'String','2');
      set( findobj(gcf,'tag','tagLower'), 'String','-2');
      set( findobj(gcf,'tag','tagFun'), 'String', ...
           'sign(x) .* sqrt( abs(x) ) + 0.5');
      set( findobj(gcf,'tag','tagDrv'),...
           'String','1 ./ ( 2*sqrt(abs(x)) )');
      set( findobj(gcf,'tag','tagEvalat'), 'String','-1');

      numnewton(4);
      warning(s);

   elseif option==50, %Set convergence options

      ConvData = get( findobj('tag','tagNewton'), 'UserData');
      prompts = {'Convergence when abs[ x_n - x_{n-1} ] is less than:',...
                 'Convergence when abs[ f(x) ] is less than:',...
                 'Maximum iterations:'};
      def = ConvData;
      dlgTitle = 'Convergence options';
      lines = 1;
      answer = inputdlg( prompts, dlgTitle, lines, def );

      %Now deal with the replies:
      if ( ~isempty(answer) ),
         set( findobj('tag','tagNewton'), 'UserData',answer);
      end

   elseif option==51, %Reset default convergence options

      ConvData = { '0.0001', '0.00001', '10' } ;
      set( findobj('tag','tagNewton'), 'UserData',ConvData);

   end;


end;

return;