function area=numdiff(option)
%NUMDIFF  Shows the numerical differentiation between points for a function
%USE:    Type  numdiff  at the MATLAB prompt

%Peter Dunn
%14 April 1999

if nargin==0,

   figure('Name','Numerical Differentiation','NumberTitle','off')

   enterframe = uicontrol('Style','Frame','Units','normalized',...
            'Position',[0.04 0.04 0.77, 0.22]);

   funtext = uicontrol('Style','text','Units','normalized',...
             'Position',[0.05 0.16 0.15,0.08],'String','Function:',...
             'FontSize',14,...
             'HorizontalAlignment','center');
   fun = uicontrol('Style','edit','Units','normalized',...
         'Position',[0.2 0.16 0.5 0.08],...
         'FontSize',14,...
         'ToolTip','The function to differentiate',...
         'String','x.^2','tag','tagDrvFun' );

   lowertext = uicontrol('Style','text','Units','normalized',...
  'FontSize',14,...
         'Position',[0.05 0.05 0.15 0.08],'String','Min x:' );
   lower = uicontrol('Style','edit','Units','normalized',...
         'FontSize',14,...
         'ToolTip','Minimum x-value to display',...
         'Position',[0.20 0.05 0.08 0.08],'String','0','tag','tagDrvLower');

   uppertext = uicontrol('Style','text','Units','normalized',...
         'FontSize',14,...
         'Position',[0.28 0.05 0.15 0.08],'String','Max x:');
   upper = uicontrol('Style','edit','Units','normalized',...
         'FontSize',14,...
         'ToolTip','Maximum x-value to display',...
         'Position',[0.43 0.05 0.08 0.08],'String','2','tag','tagDrvUpper');

   evalattext = uicontrol('Style','text','Units','normalized',...
         'FontSize',14,...
         'Position',[0.52 0.05 0.10 0.08],'String','At x=');
   evalat = uicontrol('Style','edit','Units','normalized',...
         'FontSize',14,...
         'ToolTip','Where the derivative is evaluated',...
         'Position',[0.63 0.05 0.07 0.08],'String','1','tag','tagDrvEvalat');

   speedtextup = uicontrol('Style','text','Units','normalized',...
         'FontSize',14,...
         'Position',[0.7 0.21 0.1 0.04],'String',{'Slow'});
   speedtextlo = uicontrol('Style','text','Units','normalized',...
         'FontSize',14,...
         'Position',[0.7 0.05 0.1 0.04],'String',{'Fast'});
   speed = uicontrol('Style','slider','Units','normalized',...
         'ToolTip','Changes animation speed',...
         'Position',[0.73 0.10 0.05 0.10],'Min',0,'Max',2,...
         'SliderStep',[0.2 0.50],  'Value',1,'tag','tagDrvSpeed');

   gobtn = uicontrol('Style','pushbutton','Units','normalized',...
             'Position',[0.82 0.1866 0.15 0.0733],...
             'FontSize',14,...
             'ToolTip','Start animation',...
             'String','GO','Callback','numdiff(4); ');
   pausebtn = uicontrol('Style','pushbutton','Units','normalized',...
             'Position',[0.82 0.1133 0.15 0.0733],...
             'FontSize',14,...
             'ToolTip','Pauses animation',...
             'String','PAUSE','Callback','numdiff(5);','tag','tagDrvPause');
   quitbtn = uicontrol('Style','pushbutton','Units','normalized',...
             'Position',[0.82 0.04 0.15 0.0733],...
             'FontSize',14,...
             'ToolTip','Closes program',...
             'String','QUIT','Callback','delete(gcf);');

   uimenu('Label','DEMONSTRATIONS','tag','tagDrvDemo')
   uimenu( findobj('tag', 'tagDrvDemo'), ...
           'Callback','numdiff(20);',...
           'Label','Infinite Derivative');
   uimenu( findobj('tag', 'tagDrvDemo'), ...
           'Callback','numdiff(21);',...
           'Label','Infinite Derivative (2)');
   uimenu( findobj('tag', 'tagDrvDemo'), ...
           'Callback','numdiff(22);',...
           'Label','Zero derivative');
   uimenu( findobj('tag', 'tagDrvDemo'), ...
           'Callback','numdiff(23);',...
           'Label','Zero derivative (2)');
   uimenu( findobj('tag', 'tagDrvDemo'), ...
           'Callback','numdiff(24);',...
           'Label','Difficult Derivative');
   uimenu( findobj('tag', 'tagDrvDemo'), ...
           'Callback','numdiff(25);',...
           'Label','Difficult Derivative (2)');
  

	axes('Position',[0.10 0.35 0.85, 0.6]);
   axis( [ 0 2 0 4] );
   
   set(gca,'FontSize',16);

   zoom;

   %numdiff(4);

elseif nargin==1,

   if option == 3, %NUMBER REGIONS

      xhi = str2num(get( findobj(gcf,'tag','tagDrvUpper'),'String'));
      xlo = str2num(get( findobj(gcf,'tag','tagDrvLower'),'String'));
      xval = str2num(get( findobj(gcf,'tag','tagDrvEvalat'),'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, %FUNCTION

      watchon;

      delete(gca);

      numdiff(3);

      f = get( findobj(gcf,'tag','tagDrvFun'),'String');
      
      evalat = str2num(get( findobj(gcf,'tag','tagDrvEvalat'), 'String'));

      haxes = axes('Position',[0.10 0.35 0.85, 0.6]);

      xu = str2num(get( findobj(gcf,'tag','tagDrvUpper'),'String'));
      xl = str2num(get( findobj(gcf,'tag','tagDrvLower'),'String'));

      x = linspace(xl,xu,1000);

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

      if any(~isreal(ffn)) | 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;
      plot( x, ffn, 'Linewidth',3 );
      ax = axis;
      axis(ax);

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

      plot(x, fx, 'o','MarkerFaceColor','r', ...
                      'MarkerEdgeColor','r');

      deltax = [1 0.8 0.5 0.3 0.1 0.01 0.001 0.00001 0.000000001];
      numapprox = length(deltax);

%     Display
      disp(' ' );
      disp('    Delta-x         Approximate Derivative');
      disp('    ~~~~~~~         ~~~~~~~~~~~~~~~~~~~~~~');

      for i = 1:numapprox

         
%        Evaluate at x+dx
         x = evalat+deltax(i); 
         eval( [ 'fdx = ',f,';'], 'fdx=NaN;');
      
         dx = x;

%        Tangent
         m = ( fdx - fx ) / deltax(i);
         c = fx - m*evalat;
         x = linspace(xl,xu,100);
         tangentline = m*x+c;

         if any(~isreal(fx)) | any(isnan(fx)) | ...
            any(isnan(fdx)) | 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;

         hold on;


%        Remove last tangent line
         if exist( 'htan' ), 
            delete( htan ) ; 
            delete( hothermarker );
         end;

%        Plot tangent
         htan = plot( x, tangentline, 'r-', 'LineWidth',2 );
         hothermarker = plot(dx, fdx, 'o','MarkerFaceColor','r', ...
                                         'MarkerEdgeColor','r');

%        Output
         lena = length(num2str(m));
         lenb = length(num2str(deltax(i)));
         disp([blanks(10-lenb),num2str(deltax(i)),blanks(20-lena),...
                num2str(m) ]);

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

         %Legend
         legend( htan, [ 'Slope = ',num2str( m ) ], 0 );

         drawnow;
         
         pausetime = get( findobj(gcf,'tag','tagDrvSpeed'),'value');
         pause(pausetime);

      end;

      watchoff;

   elseif option==5,

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

   elseif option==6,

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

   elseif option==20, %Infinite Derivative

      set( findobj('tag','tagDrvFun'), 'String','sign(x) .* sqrt( abs(x) )');
      set( findobj('tag','tagDrvLower'), 'String','-2');
      set( findobj('tag','tagDrvUpper'), 'String','2');
      set( findobj('tag','tagDrvEvalat'), 'String','0');
      
      numdiff(4);

   elseif  option==21, %Infinite derivative (2)

      set( findobj('tag','tagDrvFun'), 'String','besselj( 0.2, x )');
      set( findobj('tag','tagDrvLower'), 'String','0');
      set( findobj('tag','tagDrvUpper'), 'String','1');
      set( findobj('tag','tagDrvEvalat'), 'String','0');

      numdiff(4);

   elseif option==22, %Zero derivative

      set( findobj('tag','tagDrvFun'), 'String','1 - x.^2');
      set( findobj('tag','tagDrvLower'), 'String','-2');
      set( findobj('tag','tagDrvUpper'), 'String','2');
      set( findobj('tag','tagDrvEvalat'), 'String','0');

      numdiff(4);

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

      set( findobj('tag','tagDrvFun'), 'String','sin( x )');
      set( findobj('tag','tagDrvLower'), 'String','0');
      set( findobj('tag','tagDrvUpper'), 'String','2*pi');
      set( findobj('tag','tagDrvEvalat'), 'String','pi/2');

      numdiff(4);

   elseif  option==24, %Difficult derivative

      set( findobj('tag','tagDrvFun'), 'String','sin( 1./ x)');
      set( findobj('tag','tagDrvLower'), 'String','0.06');
      set( findobj('tag','tagDrvUpper'), 'String','1');
      set( findobj('tag','tagDrvEvalat'), 'String','0.09');

      numdiff(4);

   elseif  option==25, %Difficult derivative (2)

      set( findobj('tag','tagDrvFun'), 'String','humps(x)');
      set( findobj('tag','tagDrvLower'), 'String','0');
      set( findobj('tag','tagDrvUpper'), 'String','1.5');
      set( findobj('tag','tagDrvEvalat'), 'String','0.3');

      numdiff(4);

   end;


end;

return;