function plotfuns(funname,a,b,tol,details)
% Help for plotfuns.m
% copyright 1998, D. Estep
%
% PURPOSE
%
% This script plots user input functions on [a,b] on one plot at points determined
% using an adaptive strategy that insures the change in each function's
% values from one sample point to the next is smaller than a tolerance.
% The program starts with a nearly uniform mesh of sample points and then
% creates an adapted mesh by successive refinement.  A different mesh is
% computed for each function.
%
% FUNCTION STATEMENT
%
% plotfuns(funname,a,b,tol,details)
%
% REQUIRED VARIABLES
%
%    funname:     contains the name of the functions to be plotted, the functions
%                 can be built-in matlab functions and/or defined as m files
%                 
%    a,b:         numbers defining the interval [a,b] on which the function is plotted
%
% OPTIONAL VARIABLES
%
%    tol:         the tolerance for the maximum possible change in the function
%                 values at neighboring sample points. 
%    details:     indicates whether the program should give information on
%                 the sample points.
%                 details = 0 means just show plot
%                 details = 1 means show plot, changes in function, and mesh sizes
%
% USAGE
%
% First define a variable that contains the names of the functions to be plotted
% using the matlab char statement.  For instructions on using char, see help char.
% The variable is defined using a matlab statement of the form
% >> variable name = char('function1(x)','function2(x)', ...)
% where function1, function1, ... are the names of the functions, which must
% be functions of x.
%
%    For example, we would define
%    >> funname = char('sin(x)','x*exp(x)')
%    if we wanted to plot sin(x) and x*exp(x) on the same plot.
%
% A function can also be defined in a m file.  To do this, create a m 
% file named functionname.m containing the function to be plotted
% with the form:
%
% function [y] = functionname(x)
% y = give the formula in terms of input x;
%
%    For example, we create a file called quadratic.m containing the lines
%    function [y] = quadratic(x)
%    y = x^2;
%
% Now when defining the variable containing the function name, we would use
% 'functionname(x)' as one of the entries.
%
%    For the example above, we would define 
%    >> funname = char('sin(x),'quadratic(x)')
%    to get a plot of the matlab function sin(x) and the user defined
%    function quadratic(x) on the same plot.
%
% Finally, define the values of a, b, tol, details, then call
% >> plotfuns(variable containing function names, a, b, tol, details)
%
%    In the example above, we would type
%    >> plotfuns(funname,-4,4,.01,0)
%    to get a plot on the interval [-4,4] using a mesh of sample points
%    in which the change in the function values is smaller than .01
%
% The last two arguments to plotfun can be left off, for example
% >> plotfuns(funname,-4,4) 
%
%
% DETAILS
%
% The program uses two meshes: a current mesh and a new mesh.  Going from
% left to right, the difference in function values at successive nodes in
% the current mesh is checked. If this difference is smaller than tol, the
% new mesh points are defined to be the old mesh points.  If this difference
% is larger than tol, then the interval between the nodes is divided into
% uniform pieces and the new mesh is defined by using the nodes of the
% pieces.  The number of pieces is computed by assuming the function is 
% linear between the nodes of the current mesh.  After every interval in
% the current mesh is checked, the current mesh is reset to be the just
% computed new mesh.  If no refinements were needed, the iteration stops
% and if there was a refinement in some interval, there is another loop.

%initialize

% define the maximum number of adaptive iterations for any one function
maxrefineiter=5;

% check the number of arguments and set defaults for the missing values
noargs = nargin;
if noargs == 3
   tol = .1 *(b-a);
   details = 0;
elseif noargs == 4 
   details = 0;
elseif  noargs < 3  &  noargs>5 
   disp('wrong number of inputs!')
   break
end

% nf is the number of functions to be plotted
[nf lf]  = size(funname);

% loop over the functions
for jf = 1: nf
   fprintf('computing plot points for function number %d\n',jf)
   
   % define initial number of sample points and uniform sample interval
   n(jf) = 11;
   delta = (b-a)/(n(jf)-1);
   
   % define the function to be plotted, the function name is in the jf row of funname
   onefunname = funname(jf,1:lf);
   f = inline(onefunname,'x');
   
   
   % define the sample points. The first and last points align with a and b.
   x(1,jf) = a;
   y(1,jf) = f(x(1,jf));
   x(n,jf) = b;
   y(n,jf)=f(x(n(jf),jf));
   for i = 2: n(jf)-1
      x(i,jf) = a+ delta*(i-1);
      % vary the internal sample points by 10% randomly to avoid
      % problems with sampling at symmetric points
      x(i,jf) = x(i,jf) + (rand-.5)*.2*delta;
      y(i,jf)=f(x(i,jf));
   end
   
   % set refinement level
   refine = -1;
   
   % loopflag=1 means do a refinement search, the first time we do one of course
   loopflag=1;
   
   while loopflag>0 & refine <= maxrefineiter
      
      % xnew, ynew are the new mesh points and function values
      % nnew is the number of mesh points in the new mesh. These
      % are temporary variables and are not saved.
      
      % the first nodes in two meshes always agree
      xnew(1) = x(1,jf);
      ynew(1) = y(1,jf);
      nnew=1;
      
      % default: no more refinement iterations
      loopflag = -1;
      refine = refine+1;
      
      for i = 2:n(jf)
         
         % check if function difference is small enough
         if (abs(y(i,jf)-y(i-1,jf))< tol)
            
            % if so, then set new mesh point = old mesh point and go to next
            nnew=nnew+1;
            xnew(nnew) = x(i,jf);
            ynew(nnew) = y(i,jf);
            
         else
            
            % indicate that a refinement is needed in this iteration, so check
            % for one more iteration to be sure that the differences are small
            loopflag = 1;
            
            % assume that the function is linear between the two mesh points and
            % compute the necessary number of divisions to get the desired change
            nmid=ceil(abs(y(i,jf)-y(i-1,jf))/tol)+1;
            delta = abs(x(i,jf)-x(i-1,jf))/nmid;
            
            % now define the new mesh points up to the next to last one
            for k = 1:nmid-1;
               xnew(nnew+k) = x(i-1,jf) + delta*k;
               ynew(nnew+k) = f(xnew(nnew+k));
            end
            
            % treat the last point to be sure that it lines up with old mesh
            xnew(nnew+nmid) = x(i,jf);
            ynew(nnew+nmid) = y(i,jf);
            nnew = nnew + nmid;
            
         end   
         
      end
      
      % prepare for a new loop, set the old values to be the new values
      % this works because xnew and ynew are always longer or have the
      % same length as x and y.
      for i = 1: nnew
         x(i,jf) = xnew(i);
         y(i,jf) = ynew(i);
      end
      n(jf)=nnew;
      clear xnew;
      clear ynew;
      clear nnew;
      
   end
   fprintf('plotting took %d refinement steps\n', refine)

   if refine > maxrefineiter
      fprintf('The program exceeded %d refinement iterations plotting function number %d \n',maxrefineiter,jf)
      return
   end

end

% Because the function use different numbers of plotting points, some of
% the entries at the ends of the columns of x and y may be undefined.
% Fill in these entries with the last value that is defined.
M=norm(n,inf);
for jf = 1: nf
   for i = n(jf):M
      x(i,jf)=x(n(jf),jf);
      y(i,jf)=y(n(jf),jf);
   end
end

if(details==0)
   plot(x,y)
   %graphof=strcat('Graph of ',funname);
   %title('Plots of Functions')
   xlabel('x')
else
   figure
   subplot(1,3,1)
   hold on
	plot(x,y)
   %plot(x,0,'r.')
   %graphof=['Graph of ',funname];
   title('Plots')
   xlabel('x')
   hold off
   for jf = 1: nf
      for i = 1: M - 1
         dy(i,jf) = abs(y(i+1,jf)-y(i,jf));
         dx(i,jf) = abs(x(i+1,jf)-x(i,jf));
         xx(i,jf) = (x(i+1,jf)+x(i,jf))/2;
      end
   end
	subplot(1,3,2)
   plot(xx,dy)
   title('Changes in Values')
   xlabel('x')
	subplot(1,3,3)
   semilogy(xx,dx)
   title('Interval Sizes')
   xlabel('x')
	
end